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_opa_interp.F90 in branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 8902

Last change on this file since 8902 was 8902, checked in by timgraham, 6 years ago

Bug fixes for non key_vertical case, reverted to interpolating tracer rather than tracer content

  • Property svn:keywords set to Id
File size: 55.8 KB
Line 
1MODULE agrif_opa_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_opa_interp  ***
4   !! AGRIF: interpolation package
5   !!======================================================================
6   !! History :  2.0  !  2002-06  (XXX)  Original cade
7   !!             -   !  2005-11  (XXX)
8   !!            3.2  !  2009-04  (R. Benshila)
9   !!            3.6  !  2014-09  (R. Benshila)
10   !!----------------------------------------------------------------------
11#if defined key_agrif
12   !!----------------------------------------------------------------------
13   !!   'key_agrif'                                              AGRIF zoom
14   !!----------------------------------------------------------------------
15   !!   Agrif_tra     :
16   !!   Agrif_dyn     :
17   !!   interpu       :
18   !!   interpv       :
19   !!----------------------------------------------------------------------
20   USE par_oce
21   USE oce
22   USE dom_oce     
23   USE zdf_oce
24   USE agrif_oce
25   USE phycst
26   !
27   USE in_out_manager
28   USE agrif_opa_sponge
29   USE lib_mpp
30   USE wrk_nemo
31 
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
36   PUBLIC   interpun, interpvn
37   PUBLIC   interptsn,  interpsshn
38   PUBLIC   interpunb, interpvnb, interpub2b, interpvb2b
39   PUBLIC   interpe3t, interpumsk, interpvmsk
40# if defined key_zdftke
41   PUBLIC   Agrif_tke, interpavm
42# endif
43
44   INTEGER ::   bdy_tinterp = 0
45
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !! NEMO/NST 3.7 , NEMO Consortium (2015)
49   !! $Id$
50   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE Agrif_tra
55      !!----------------------------------------------------------------------
56      !!                  ***  ROUTINE Agrif_tra  ***
57      !!----------------------------------------------------------------------
58      !
59      IF( Agrif_Root() )   RETURN
60      !
61      Agrif_SpecialValue    = 0._wp
62      Agrif_UseSpecialValue = .TRUE.
63      !
64      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
65      !
66      Agrif_UseSpecialValue = .FALSE.
67      !
68   END SUBROUTINE Agrif_tra
69
70
71   SUBROUTINE Agrif_dyn( kt )
72      !!----------------------------------------------------------------------
73      !!                  ***  ROUTINE Agrif_DYN  ***
74      !!---------------------------------------------------------------------- 
75      INTEGER, INTENT(in) ::   kt
76      !
77      INTEGER ::   ji, jj, jk       ! dummy loop indices
78      INTEGER ::   j1, j2, i1, i2
79      REAL(wp), POINTER, DIMENSION(:,:) ::   zub, zvb
80      !!---------------------------------------------------------------------- 
81      !
82      IF( Agrif_Root() )   RETURN
83      !
84      CALL wrk_alloc( jpi,jpj,   zub, zvb )
85      !
86      Agrif_SpecialValue    = 0._wp
87      Agrif_UseSpecialValue = ln_spc_dyn
88      !
89      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
90      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
91      !
92      Agrif_UseSpecialValue = .FALSE.
93      !
94      ! prevent smoothing in ghost cells
95      i1 =  1   ;   i2 = jpi
96      j1 =  1   ;   j2 = jpj
97      IF( nbondj == -1 .OR. nbondj == 2 )   j1 = 3
98      IF( nbondj == +1 .OR. nbondj == 2 )   j2 = nlcj-2
99      IF( nbondi == -1 .OR. nbondi == 2 )   i1 = 3
100      IF( nbondi == +1 .OR. nbondi == 2 )   i2 = nlci-2
101
102      IF( nbondi == -1 .OR. nbondi == 2 ) THEN
103         !
104         ! Smoothing
105         ! ---------
106         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
107            ua_b(2,:) = 0._wp
108            DO jk = 1, jpkm1
109               DO jj = 1, jpj
110                  ua_b(2,jj) = ua_b(2,jj) + e3u_a(2,jj,jk) * ua(2,jj,jk)
111               END DO
112            END DO
113            DO jj = 1, jpj
114               ua_b(2,jj) = ua_b(2,jj) * r1_hu_a(2,jj)           
115            END DO
116         ENDIF
117         !
118         DO jk=1,jpkm1                 ! Smooth
119            DO jj=j1,j2
120               ua(2,jj,jk) = 0.25_wp*(ua(1,jj,jk)+2._wp*ua(2,jj,jk)+ua(3,jj,jk))
121               ua(2,jj,jk) = ua(2,jj,jk) * umask(2,jj,jk)
122            END DO
123         END DO
124         !
125         zub(2,:) = 0._wp              ! Correct transport
126         DO jk = 1, jpkm1
127            DO jj = 1, jpj
128               zub(2,jj) = zub(2,jj) + e3u_a(2,jj,jk) * ua(2,jj,jk)
129            END DO
130         END DO
131         DO jj=1,jpj
132            zub(2,jj) = zub(2,jj) * r1_hu_a(2,jj)
133         END DO
134
135         DO jk=1,jpkm1
136            DO jj=1,jpj
137               ua(2,jj,jk) = (ua(2,jj,jk)+ua_b(2,jj)-zub(2,jj))*umask(2,jj,jk)
138            END DO
139         END DO
140
141         ! Set tangential velocities to time splitting estimate
142         !-----------------------------------------------------
143         IF( ln_dynspg_ts ) THEN
144            zvb(2,:) = 0._wp
145            DO jk = 1, jpkm1
146               DO jj = 1, jpj
147                  zvb(2,jj) = zvb(2,jj) + e3v_a(2,jj,jk) * va(2,jj,jk)
148               END DO
149            END DO
150            DO jj = 1, jpj
151               zvb(2,jj) = zvb(2,jj) * r1_hv_a(2,jj)
152            END DO
153            DO jk = 1, jpkm1
154               DO jj = 1, jpj
155                  va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-zvb(2,jj)) * vmask(2,jj,jk)
156               END DO
157            END DO
158         ENDIF
159         !
160         ! Mask domain edges:
161         !-------------------
162         DO jk = 1, jpkm1
163            DO jj = 1, jpj
164               ua(1,jj,jk) = 0._wp
165               va(1,jj,jk) = 0._wp
166            END DO
167         END DO         
168         !
169      ENDIF
170
171      IF( nbondi == 1 .OR. nbondi == 2 ) THEN
172
173         ! Smoothing
174         ! ---------
175         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
176            ua_b(nlci-2,:) = 0._wp
177            DO jk=1,jpkm1
178               DO jj=1,jpj
179                  ua_b(nlci-2,jj) = ua_b(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
180               END DO
181            END DO
182            DO jj=1,jpj
183               ua_b(nlci-2,jj) = ua_b(nlci-2,jj) * r1_hu_a(nlci-2,jj)           
184            END DO
185         ENDIF
186
187         DO jk = 1, jpkm1              ! Smooth
188            DO jj = j1, j2
189               ua(nlci-2,jj,jk) = 0.25_wp * umask(nlci-2,jj,jk)      &
190                  &             * ( ua(nlci-3,jj,jk) + 2._wp*ua(nlci-2,jj,jk) + ua(nlci-1,jj,jk) )
191            END DO
192         END DO
193
194         zub(nlci-2,:) = 0._wp        ! Correct transport
195         DO jk = 1, jpkm1
196            DO jj = 1, jpj
197               zub(nlci-2,jj) = zub(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
198            END DO
199         END DO
200         DO jj = 1, jpj
201            zub(nlci-2,jj) = zub(nlci-2,jj) * r1_hu_a(nlci-2,jj)
202         END DO
203
204         DO jk = 1, jpkm1
205            DO jj = 1, jpj
206               ua(nlci-2,jj,jk) = ( ua(nlci-2,jj,jk) + ua_b(nlci-2,jj) - zub(nlci-2,jj) ) * umask(nlci-2,jj,jk)
207            END DO
208         END DO
209         !
210         ! Set tangential velocities to time splitting estimate
211         !-----------------------------------------------------
212         IF( ln_dynspg_ts ) THEN
213            zvb(nlci-1,:) = 0._wp
214            DO jk = 1, jpkm1
215               DO jj = 1, jpj
216                  zvb(nlci-1,jj) = zvb(nlci-1,jj) + e3v_a(nlci-1,jj,jk) * va(nlci-1,jj,jk)
217               END DO
218            END DO
219            DO jj=1,jpj
220               zvb(nlci-1,jj) = zvb(nlci-1,jj) * r1_hv_a(nlci-1,jj)
221            END DO
222            DO jk = 1, jpkm1
223               DO jj = 1, jpj
224                  va(nlci-1,jj,jk) = ( va(nlci-1,jj,jk) + va_b(nlci-1,jj) - zvb(nlci-1,jj) ) * vmask(nlci-1,jj,jk)
225               END DO
226            END DO
227         ENDIF
228         !
229         ! Mask domain edges:
230         !-------------------
231         DO jk = 1, jpkm1
232            DO jj = 1, jpj
233               ua(nlci-1,jj,jk) = 0._wp
234               va(nlci  ,jj,jk) = 0._wp
235            END DO
236         END DO 
237         !
238      ENDIF
239
240      IF( nbondj == -1 .OR. nbondj == 2 ) THEN
241
242         ! Smoothing
243         ! ---------
244         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
245            va_b(:,2) = 0._wp
246            DO jk = 1, jpkm1
247               DO ji = 1, jpi
248                  va_b(ji,2) = va_b(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk)
249               END DO
250            END DO
251            DO ji=1,jpi
252               va_b(ji,2) = va_b(ji,2) * r1_hv_a(ji,2)           
253            END DO
254         ENDIF
255         !
256         DO jk = 1, jpkm1              ! Smooth
257            DO ji = i1, i2
258               va(ji,2,jk) = 0.25_wp * vmask(ji,2,jk)    &
259                  &        * ( va(ji,1,jk) + 2._wp*va(ji,2,jk) + va(ji,3,jk) )
260            END DO
261         END DO
262         !
263         zvb(:,2) = 0._wp              ! Correct transport
264         DO jk=1,jpkm1
265            DO ji=1,jpi
266               zvb(ji,2) = zvb(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk)
267            END DO
268         END DO
269         DO ji = 1, jpi
270            zvb(ji,2) = zvb(ji,2) * r1_hv_a(ji,2)
271         END DO
272         DO jk = 1, jpkm1
273            DO ji = 1, jpi
274               va(ji,2,jk) = ( va(ji,2,jk) + va_b(ji,2) - zvb(ji,2) ) * vmask(ji,2,jk)
275            END DO
276         END DO
277
278         ! Set tangential velocities to time splitting estimate
279         !-----------------------------------------------------
280         IF( ln_dynspg_ts ) THEN
281            zub(:,2) = 0._wp
282            DO jk = 1, jpkm1
283               DO ji = 1, jpi
284                  zub(ji,2) = zub(ji,2) + e3u_a(ji,2,jk) * ua(ji,2,jk) * umask(ji,2,jk)
285               END DO
286            END DO
287            DO ji = 1, jpi
288               zub(ji,2) = zub(ji,2) * r1_hu_a(ji,2)
289            END DO
290
291            DO jk = 1, jpkm1
292               DO ji = 1, jpi
293                  ua(ji,2,jk) = ( ua(ji,2,jk) + ua_b(ji,2) - zub(ji,2) ) * umask(ji,2,jk)
294               END DO
295            END DO
296         ENDIF
297
298         ! Mask domain edges:
299         !-------------------
300         DO jk = 1, jpkm1
301            DO ji = 1, jpi
302               ua(ji,1,jk) = 0._wp
303               va(ji,1,jk) = 0._wp
304            END DO
305         END DO
306
307      ENDIF
308
309      IF( nbondj == 1 .OR. nbondj == 2 ) THEN
310         !
311         ! Smoothing
312         ! ---------
313         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
314            va_b(:,nlcj-2) = 0._wp
315            DO jk = 1, jpkm1
316               DO ji = 1, jpi
317                  va_b(ji,nlcj-2) = va_b(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk)
318               END DO
319            END DO
320            DO ji = 1, jpi
321               va_b(ji,nlcj-2) = va_b(ji,nlcj-2) * r1_hv_a(ji,nlcj-2)           
322            END DO
323         ENDIF
324         !
325         DO jk = 1, jpkm1              ! Smooth
326            DO ji = i1, i2
327               va(ji,nlcj-2,jk) = 0.25_wp * vmask(ji,nlcj-2,jk)   &
328                  &             * ( va(ji,nlcj-3,jk) + 2._wp * va(ji,nlcj-2,jk) + va(ji,nlcj-1,jk) )
329            END DO
330         END DO
331         !
332         zvb(:,nlcj-2) = 0._wp         ! Correct transport
333         DO jk = 1, jpkm1
334            DO ji = 1, jpi
335               zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk)
336            END DO
337         END DO
338         DO ji = 1, jpi
339            zvb(ji,nlcj-2) = zvb(ji,nlcj-2) * r1_hv_a(ji,nlcj-2)
340         END DO
341         DO jk = 1, jpkm1
342            DO ji = 1, jpi
343               va(ji,nlcj-2,jk) = ( va(ji,nlcj-2,jk) + va_b(ji,nlcj-2) - zvb(ji,nlcj-2) ) * vmask(ji,nlcj-2,jk)
344            END DO
345         END DO
346         !
347         ! Set tangential velocities to time splitting estimate
348         !-----------------------------------------------------
349         IF( ln_dynspg_ts ) THEN
350            zub(:,nlcj-1) = 0._wp
351            DO jk = 1, jpkm1
352               DO ji = 1, jpi
353                  zub(ji,nlcj-1) = zub(ji,nlcj-1) + e3u_a(ji,nlcj-1,jk) * ua(ji,nlcj-1,jk) * umask(ji,nlcj-1,jk)
354               END DO
355            END DO
356            DO ji = 1, jpi
357               zub(ji,nlcj-1) = zub(ji,nlcj-1) * r1_hu_a(ji,nlcj-1)
358            END DO
359            !
360            DO jk = 1, jpkm1
361               DO ji = 1, jpi
362                  ua(ji,nlcj-1,jk) = ( ua(ji,nlcj-1,jk) + ua_b(ji,nlcj-1) - zub(ji,nlcj-1) ) * umask(ji,nlcj-1,jk)
363               END DO
364            END DO
365         ENDIF
366         !
367         ! Mask domain edges:
368         !-------------------
369         DO jk = 1, jpkm1
370            DO ji = 1, jpi
371               ua(ji,nlcj  ,jk) = 0._wp
372               va(ji,nlcj-1,jk) = 0._wp
373            END DO
374         END DO 
375         !
376      ENDIF
377      !
378      CALL wrk_dealloc( jpi,jpj,   zub, zvb )
379      !
380   END SUBROUTINE Agrif_dyn
381
382
383   SUBROUTINE Agrif_dyn_ts( jn )
384      !!----------------------------------------------------------------------
385      !!                  ***  ROUTINE Agrif_dyn_ts  ***
386      !!---------------------------------------------------------------------- 
387      !!
388      INTEGER, INTENT(in) ::   jn
389      !!
390      INTEGER :: ji, jj
391      !!---------------------------------------------------------------------- 
392      !
393      IF( Agrif_Root() )   RETURN
394      !
395      IF((nbondi == -1).OR.(nbondi == 2)) THEN
396         DO jj=1,jpj
397            va_e(2,jj) = vbdy_w(jj) * hvr_e(2,jj)
398            ! Specified fluxes:
399            ua_e(2,jj) = ubdy_w(jj) * hur_e(2,jj)
400            ! Characteristics method:
401            !alt            ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) &
402            !alt                       &           - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) )
403         END DO
404      ENDIF
405      !
406      IF((nbondi == 1).OR.(nbondi == 2)) THEN
407         DO jj=1,jpj
408            va_e(nlci-1,jj) = vbdy_e(jj) * hvr_e(nlci-1,jj)
409            ! Specified fluxes:
410            ua_e(nlci-2,jj) = ubdy_e(jj) * hur_e(nlci-2,jj)
411            ! Characteristics method:
412            !alt            ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) &
413            !alt                            &           + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) )
414         END DO
415      ENDIF
416      !
417      IF((nbondj == -1).OR.(nbondj == 2)) THEN
418         DO ji=1,jpi
419            ua_e(ji,2) = ubdy_s(ji) * hur_e(ji,2)
420            ! Specified fluxes:
421            va_e(ji,2) = vbdy_s(ji) * hvr_e(ji,2)
422            ! Characteristics method:
423            !alt            va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) &
424            !alt                       &           - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) )
425         END DO
426      ENDIF
427      !
428      IF((nbondj == 1).OR.(nbondj == 2)) THEN
429         DO ji=1,jpi
430            ua_e(ji,nlcj-1) = ubdy_n(ji) * hur_e(ji,nlcj-1)
431            ! Specified fluxes:
432            va_e(ji,nlcj-2) = vbdy_n(ji) * hvr_e(ji,nlcj-2)
433            ! Characteristics method:
434            !alt            va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2)  + va_e(ji,nlcj-3) &
435            !alt                            &           + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) )
436         END DO
437      ENDIF
438      !
439   END SUBROUTINE Agrif_dyn_ts
440
441
442   SUBROUTINE Agrif_dta_ts( kt )
443      !!----------------------------------------------------------------------
444      !!                  ***  ROUTINE Agrif_dta_ts  ***
445      !!---------------------------------------------------------------------- 
446      !!
447      INTEGER, INTENT(in) ::   kt
448      !!
449      INTEGER :: ji, jj
450      LOGICAL :: ll_int_cons
451      REAL(wp) :: zrhot, zt
452      !!---------------------------------------------------------------------- 
453      !
454      IF( Agrif_Root() )   RETURN
455      !
456      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
457      !
458      zrhot = Agrif_rhot()
459      !
460      ! "Central" time index for interpolation:
461      IF( ln_bt_fw ) THEN
462         zt = REAL( Agrif_NbStepint()+0.5_wp, wp ) / zrhot
463      ELSE
464         zt = REAL( Agrif_NbStepint()       , wp ) / zrhot
465      ENDIF
466      !
467      ! Linear interpolation of sea level
468      Agrif_SpecialValue    = 0._wp
469      Agrif_UseSpecialValue = .TRUE.
470      CALL Agrif_Bc_variable( sshn_id, calledweight=zt, procname=interpsshn )
471      Agrif_UseSpecialValue = .FALSE.
472      !
473      ! Interpolate barotropic fluxes
474      Agrif_SpecialValue=0.
475      Agrif_UseSpecialValue = ln_spc_dyn
476      !
477      IF( ll_int_cons ) THEN  ! Conservative interpolation
478         ! orders matters here !!!!!!
479         CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
480         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )
481         bdy_tinterp = 1
482         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
483         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )
484         bdy_tinterp = 2
485         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
486         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )         
487      ELSE ! Linear interpolation
488         bdy_tinterp = 0
489         ubdy_w(:) = 0._wp   ;   vbdy_w(:) = 0._wp 
490         ubdy_e(:) = 0._wp   ;   vbdy_e(:) = 0._wp 
491         ubdy_n(:) = 0._wp   ;   vbdy_n(:) = 0._wp 
492         ubdy_s(:) = 0._wp   ;   vbdy_s(:) = 0._wp
493         CALL Agrif_Bc_variable( unb_id, calledweight=zt, procname=interpunb )
494         CALL Agrif_Bc_variable( vnb_id, calledweight=zt, procname=interpvnb )
495      ENDIF
496      Agrif_UseSpecialValue = .FALSE.
497      !
498   END SUBROUTINE Agrif_dta_ts
499
500
501   SUBROUTINE Agrif_ssh( kt )
502      !!----------------------------------------------------------------------
503      !!                  ***  ROUTINE Agrif_DYN  ***
504      !!---------------------------------------------------------------------- 
505      INTEGER, INTENT(in) ::   kt
506      !!
507      !!---------------------------------------------------------------------- 
508      !
509      IF( Agrif_Root() )   RETURN
510      !
511      IF((nbondi == -1).OR.(nbondi == 2)) THEN
512         ssha(2,:)=ssha(3,:)
513         sshn(2,:)=sshn(3,:)
514      ENDIF
515      !
516      IF((nbondi == 1).OR.(nbondi == 2)) THEN
517         ssha(nlci-1,:)=ssha(nlci-2,:)
518         sshn(nlci-1,:)=sshn(nlci-2,:)
519      ENDIF
520      !
521      IF((nbondj == -1).OR.(nbondj == 2)) THEN
522         ssha(:,2)=ssha(:,3)
523         sshn(:,2)=sshn(:,3)
524      ENDIF
525      !
526      IF((nbondj == 1).OR.(nbondj == 2)) THEN
527         ssha(:,nlcj-1)=ssha(:,nlcj-2)
528         sshn(:,nlcj-1)=sshn(:,nlcj-2)
529      ENDIF
530      !
531   END SUBROUTINE Agrif_ssh
532
533
534   SUBROUTINE Agrif_ssh_ts( jn )
535      !!----------------------------------------------------------------------
536      !!                  ***  ROUTINE Agrif_ssh_ts  ***
537      !!---------------------------------------------------------------------- 
538      INTEGER, INTENT(in) ::   jn
539      !!
540      INTEGER :: ji,jj
541      !!---------------------------------------------------------------------- 
542      !
543      IF((nbondi == -1).OR.(nbondi == 2)) THEN
544         DO jj = 1, jpj
545            ssha_e(2,jj) = hbdy_w(jj)
546         END DO
547      ENDIF
548      !
549      IF((nbondi == 1).OR.(nbondi == 2)) THEN
550         DO jj = 1, jpj
551            ssha_e(nlci-1,jj) = hbdy_e(jj)
552         END DO
553      ENDIF
554      !
555      IF((nbondj == -1).OR.(nbondj == 2)) THEN
556         DO ji = 1, jpi
557            ssha_e(ji,2) = hbdy_s(ji)
558         END DO
559      ENDIF
560      !
561      IF((nbondj == 1).OR.(nbondj == 2)) THEN
562         DO ji = 1, jpi
563            ssha_e(ji,nlcj-1) = hbdy_n(ji)
564         END DO
565      ENDIF
566      !
567   END SUBROUTINE Agrif_ssh_ts
568
569# if defined key_zdftke
570
571   SUBROUTINE Agrif_tke
572      !!----------------------------------------------------------------------
573      !!                  ***  ROUTINE Agrif_tke  ***
574      !!---------------------------------------------------------------------- 
575      REAL(wp) ::   zalpha
576      !!---------------------------------------------------------------------- 
577      !
578     
579      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
580      IF( zalpha > 1. )   zalpha = 1.
581      !
582      Agrif_SpecialValue    = 0.e0
583      Agrif_UseSpecialValue = .TRUE.
584      !
585      CALL Agrif_Bc_variable(avm_id ,calledweight=zalpha, procname=interpavm)       
586      !
587      Agrif_UseSpecialValue = .FALSE.
588      !
589   END SUBROUTINE Agrif_tke
590   
591# endif
592
593   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
594      !!----------------------------------------------------------------------
595      !!   *** ROUTINE interptsn ***
596      !!----------------------------------------------------------------------
597      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
598      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
599      LOGICAL                                     , INTENT(in   ) ::   before
600      INTEGER                                     , INTENT(in   ) ::   nb , ndir
601      !
602      INTEGER  ::   ji, jj, jk, jn, iref, jref   ! dummy loop indices
603      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
604      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
605      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
606      LOGICAL :: western_side, eastern_side,northern_side,southern_side
607      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
608      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
609      REAL(wp) :: h_in(k1:k2)
610      REAL(wp) :: h_out(1:jpk)
611      REAL(wp) :: h_diff, zrhoxy
612
613      zrhoxy = Agrif_rhox()*Agrif_rhoy()
614      IF (before) THEN         
615         DO jn = n1,n2-1
616            DO jk=k1,k2
617               DO jj=j1,j2
618                 DO ji=i1,i2
619                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
620                 END DO
621              END DO
622           END DO
623        END DO
624        DO jk=k1,k2
625           DO jj=j1,j2
626              DO ji=i1,i2
627                 ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
628              END DO
629           END DO
630        END DO
631      ELSE
632         western_side  = (nb == 1).AND.(ndir == 1)
633         eastern_side  = (nb == 1).AND.(ndir == 2)
634         southern_side = (nb == 2).AND.(ndir == 1)
635         northern_side = (nb == 2).AND.(ndir == 2)
636#ifdef key_vertical         
637         ptab_child(:,:,:,:) = 0.
638         DO jj=j1,j2
639            DO ji=i1,i2
640               iref = ji
641               jref = jj
642               if(western_side) iref=2
643               if(eastern_side) iref=nlci-1
644               if(southern_side) jref=2
645               if(northern_side) jref=nlcj-1
646               N_in = 0
647               DO jk=k1,k2 !k2 = jpk of parent grid
648                  IF (ptab(ji,jj,jk,n2) == 0) EXIT
649                  N_in = N_in + 1
650                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
651                  h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy)
652               END DO
653               N_out = 0
654               DO jk=1,jpk ! jpk of child grid
655                  IF (tmask(iref,jref,jk) == 0) EXIT
656                  N_out = N_out + 1
657                  h_out(jk) = e3t_n(iref,jref,jk)
658               ENDDO
659               IF (N_in > 0) THEN
660                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
661                  DO jn=1,jpts
662                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
663                  ENDDO
664               ENDIF
665            ENDDO
666         ENDDO
667#else
668         DO jk=k1,k2
669            DO jj=j1,j2
670               DO ji=i1,i2
671                  ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)
672               ENDDO
673            ENDDO
674         ENDDO
675#endif
676         !
677         zrhox = Agrif_Rhox()
678         !
679         zalpha1 = ( zrhox - 1. ) * 0.5
680         zalpha2 = 1. - zalpha1
681         !
682         zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. )
683         zalpha4 = 1. - zalpha3
684         !
685         zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
686         zalpha7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
687         zalpha5 = 1. - zalpha6 - zalpha7
688         !
689         imin = i1
690         imax = i2
691         jmin = j1
692         jmax = j2
693         !
694         ! Remove CORNERS
695         IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3
696         IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2
697         IF((nbondi == -1).OR.(nbondi == 2)) imin = 3
698         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2       
699         !
700         IF( eastern_side ) THEN
701            DO jn = 1, jpts
702               tsa(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn)
703               DO jk = 1, jpkm1
704                  DO jj = jmin,jmax
705                     IF( umask(nlci-2,jj,jk) == 0._wp ) THEN
706                        tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
707                     ELSE
708                        tsa(nlci-1,jj,jk,jn)=(zalpha4*tsa(nlci,jj,jk,jn)+zalpha3*tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
709                        IF( un(nlci-2,jj,jk) > 0._wp ) THEN
710                           tsa(nlci-1,jj,jk,jn)=( zalpha6*tsa(nlci-2,jj,jk,jn)+zalpha5*tsa(nlci,jj,jk,jn) & 
711                                 + zalpha7*tsa(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk)
712                        ENDIF
713                     ENDIF
714                  END DO
715               END DO
716               tsa(nlci,j1:j2,k1:k2,jn) = 0._wp
717            END DO
718         ENDIF
719         !
720         IF( northern_side ) THEN           
721            DO jn = 1, jpts
722               tsa(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn)
723               DO jk = 1, jpkm1
724                  DO ji = imin,imax
725                     IF( vmask(ji,nlcj-2,jk) == 0._wp ) THEN
726                        tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
727                     ELSE
728                        tsa(ji,nlcj-1,jk,jn)=(zalpha4*tsa(ji,nlcj,jk,jn)+zalpha3*tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)       
729                        IF (vn(ji,nlcj-2,jk) > 0._wp ) THEN
730                           tsa(ji,nlcj-1,jk,jn)=( zalpha6*tsa(ji,nlcj-2,jk,jn)+zalpha5*tsa(ji,nlcj,jk,jn)  &
731                                 + zalpha7*tsa(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk)
732                        ENDIF
733                     ENDIF
734                  END DO
735               END DO
736               tsa(i1:i2,nlcj,k1:k2,jn) = 0._wp
737            END DO
738         ENDIF
739         !
740         IF( western_side ) THEN           
741            DO jn = 1, jpts
742               tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn)
743               DO jk = 1, jpkm1
744                  DO jj = jmin,jmax
745                     IF( umask(2,jj,jk) == 0._wp ) THEN
746                        tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
747                     ELSE
748                        tsa(2,jj,jk,jn)=(zalpha4*tsa(1,jj,jk,jn)+zalpha3*tsa(3,jj,jk,jn))*tmask(2,jj,jk)       
749                        IF( un(2,jj,jk) < 0._wp ) THEN
750                           tsa(2,jj,jk,jn)=(zalpha6*tsa(3,jj,jk,jn)+zalpha5*tsa(1,jj,jk,jn)+zalpha7*tsa(4,jj,jk,jn))*tmask(2,jj,jk)
751                        ENDIF
752                     ENDIF
753                  END DO
754               END DO
755               tsa(1,j1:j2,k1:k2,jn) = 0._wp
756            END DO
757         ENDIF
758         !
759         IF( southern_side ) THEN           
760            DO jn = 1, jpts
761               tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn)
762               DO jk = 1, jpk     
763                  DO ji=imin,imax
764                     IF( vmask(ji,2,jk) == 0._wp ) THEN
765                        tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
766                     ELSE
767                        tsa(ji,2,jk,jn)=(zalpha4*tsa(ji,1,jk,jn)+zalpha3*tsa(ji,3,jk,jn))*tmask(ji,2,jk)
768                        IF( vn(ji,2,jk) < 0._wp ) THEN
769                           tsa(ji,2,jk,jn)=(zalpha6*tsa(ji,3,jk,jn)+zalpha5*tsa(ji,1,jk,jn)+zalpha7*tsa(ji,4,jk,jn))*tmask(ji,2,jk)
770                        ENDIF
771                     ENDIF
772                  END DO
773               END DO
774               tsa(i1:i2,1,k1:k2,jn) = 0._wp
775            END DO
776         ENDIF
777         !
778         ! Treatment of corners
779         !
780         ! East south
781         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
782            tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts)
783         ENDIF
784         ! East north
785         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
786            tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts)
787         ENDIF
788         ! West south
789         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
790            tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts)
791         ENDIF
792         ! West north
793         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
794            tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts)
795         ENDIF
796         !
797      ENDIF
798      !
799   END SUBROUTINE interptsn
800
801
802   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
803      !!----------------------------------------------------------------------
804      !!                  ***  ROUTINE interpsshn  ***
805      !!---------------------------------------------------------------------- 
806      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
807      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
808      LOGICAL                         , INTENT(in   ) ::   before
809      INTEGER                         , INTENT(in   ) ::   nb , ndir
810      !
811      LOGICAL :: western_side, eastern_side,northern_side,southern_side
812      !!---------------------------------------------------------------------- 
813      !
814      IF( before) THEN
815         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
816      ELSE
817         western_side  = (nb == 1).AND.(ndir == 1)
818         eastern_side  = (nb == 1).AND.(ndir == 2)
819         southern_side = (nb == 2).AND.(ndir == 1)
820         northern_side = (nb == 2).AND.(ndir == 2)
821         IF(western_side)  hbdy_w(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
822         IF(eastern_side)  hbdy_e(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
823         IF(southern_side) hbdy_s(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
824         IF(northern_side) hbdy_n(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
825      ENDIF
826      !
827   END SUBROUTINE interpsshn
828
829
830   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
831      !!----------------------------------------------------------------------
832      !!   *** ROUTINE interpun ***
833      !!---------------------------------------------   
834      !!
835      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
836      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
837      LOGICAL, INTENT(in) :: before
838      INTEGER, INTENT(in) :: nb , ndir
839      !!
840      INTEGER :: ji,jj,jk
841      REAL(wp) :: zrhoy
842! VERTICAL REFINEMENT BEGIN
843      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
844      REAL(wp), DIMENSION(k1:k2) :: tabin
845      REAL(wp) :: h_in(k1:k2)
846      REAL(wp) :: h_out(1:jpk)
847      INTEGER :: N_in, N_out
848      REAL(wp) :: h_diff
849      LOGICAL :: western_side, eastern_side
850      INTEGER :: iref
851
852! VERTICAL REFINEMENT END
853      !!---------------------------------------------   
854      !
855      zrhoy = Agrif_rhoy()
856      IF (before) THEN 
857         !We can't use zero as the special value because we need to include zeros
858         !when interpolating the scale factors
859         IF(Agrif_UseSpecialValue) THEN 
860!             Agrif_SpecialValue = -999._wp
861             Agrif_SpecialValue = 0._wp
862         ELSE
863             Agrif_SpecialValue = 0._wp
864         ENDIF
865         DO jk=1,jpk
866            DO jj=j1,j2
867               DO ji=i1,i2
868                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) - &
869                                   & ((umask(ji,jj,jk)-1) * Agrif_SpecialValue)
870                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
871               END DO
872            END DO
873         END DO
874      ELSE
875! VERTICAL REFINEMENT BEGIN
876         western_side  = (nb == 1).AND.(ndir == 1)
877         eastern_side  = (nb == 1).AND.(ndir == 2)
878
879         Agrif_SpecialValue = 0._wp ! reset specialvalue to zero now interpolation completed
880
881         ptab_child(:,:,:) = 0.
882#ifdef key_vertical
883! VERTICAL REFINEMENT BEGIN
884         DO jj=j1,j2
885            DO ji=i1,i2
886               iref = ji
887               IF (western_side) iref = 2
888               IF (eastern_side) iref = nlci-2
889
890               N_in = 0
891               DO jk=k1,k2
892                  IF (ptab(ji,jj,jk,2) == 0) EXIT
893                  N_in = N_in + 1
894                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
895                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
896              ENDDO
897         
898              IF (N_in == 0) THEN
899                 ptab_child(ji,jj,:) = 0.
900                 CYCLE
901              ENDIF
902         
903              N_out = 0
904              DO jk=1,jpk
905                 if (umask(ji,jj,jk) == 0) EXIT
906                 N_out = N_out + 1
907                 h_out(N_out) = e3u_n(ji,jj,jk)
908              ENDDO
909         
910              IF (N_out == 0) THEN
911                 ptab_child(ji,jj,:) = 0.
912                 CYCLE
913              ENDIF
914         
915              IF (N_in * N_out > 0) THEN
916                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
917! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
918                 if (h_diff < -1.e4) then
919                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
920!                    stop
921                 endif
922              ENDIF
923              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
924            ENDDO
925         ENDDO
926
927! VERTICAL REFINEMENT END
928         DO jk = 1, jpkm1
929            DO jj=j1,j2
930               ua(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
931            END DO
932         END DO
933#else
934         DO jk = 1, jpkm1
935            DO jj=j1,j2
936               ua(i1:i2,jj,jk) = umask(i1:i2,jj,jk) * ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_n(i1:i2,jj,jk) )
937            END DO
938         END DO
939#endif
940      ENDIF
941      !
942   END SUBROUTINE interpun
943
944
945   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
946      !!----------------------------------------------------------------------
947      !!   *** ROUTINE interpvn ***
948      !!----------------------------------------------------------------------
949      !
950      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
951      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
952      LOGICAL, INTENT(in) :: before
953      INTEGER, INTENT(in) :: nb , ndir
954      !
955      INTEGER :: ji,jj,jk
956      REAL(wp) :: zrhox
957      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
958      REAL(wp), DIMENSION(k1:k2) :: tabin
959      REAL(wp) :: h_in(k1:k2)
960      REAL(wp) :: h_out(1:jpk)
961      INTEGER :: N_in, N_out
962      REAL(wp) :: h_diff
963      LOGICAL :: northern_side,southern_side
964      INTEGER :: jref
965
966      !!---------------------------------------------   
967      !     
968      zrhox = Agrif_rhox()
969      IF (before) THEN         
970         IF(Agrif_UseSpecialValue) THEN 
971!             Agrif_SpecialValue = -999._wp
972             Agrif_SpecialValue = 0._wp
973         ELSE
974             Agrif_SpecialValue = 0._wp
975         ENDIF
976         DO jk=k1,k2
977            DO jj=j1,j2
978               DO ji=i1,i2
979                  ptab(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
980                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
981               END DO
982            END DO
983         END DO
984      ELSE       
985         Agrif_SpecialValue = 0._wp !Reset special value to zero now interpolation is done
986         ptab_child(:,:,:) = 0.
987#ifdef key_vertical
988         southern_side = (nb == 2).AND.(ndir == 1)
989         northern_side = (nb == 2).AND.(ndir == 2)
990
991         do jj=j1,j2
992            jref = jj
993            IF (southern_side) jref = 2
994            IF (northern_side) jref = nlcj-2
995            do ji=i1,i2
996
997               N_in = 0
998               do jk=k1,k2
999                  if (ptab(ji,jj,jk,2) == 0) EXIT
1000                  N_in = N_in + 1
1001                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1002                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
1003               enddo
1004               IF (N_in == 0) THEN
1005                  ptab_child(ji,jj,:) = 0.
1006                  CYCLE
1007               ENDIF
1008         
1009               N_out = 0
1010               do jk=1,jpk
1011                  if (vmask(ji,jref,jk) == 0) EXIT
1012                  N_out = N_out + 1
1013                  h_out(N_out) = e3v_n(ji,jj,jk)
1014               enddo
1015               IF (N_out == 0) THEN
1016                 ptab_child(ji,jj,:) = 0.
1017                 CYCLE
1018               ENDIF
1019               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
1020            enddo
1021         enddo
1022            DO jk=1,jpkm1
1023               DO jj=j1,j2
1024                  va(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
1025               END DO
1026            END DO
1027#else
1028            DO jk=1,jpkm1
1029               va(i1:i2,j1:j2,jk) = vmask(i1:i2,j1:j2,jk) * ptab(i1:i2,j1:j2,jk,1) / &
1030                                    & ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) )
1031            END DO
1032#endif
1033      ENDIF
1034      !       
1035   END SUBROUTINE interpvn
1036   
1037
1038   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
1039      !!----------------------------------------------------------------------
1040      !!                  ***  ROUTINE interpunb  ***
1041      !!---------------------------------------------------------------------- 
1042      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1043      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1044      LOGICAL                         , INTENT(in   ) ::   before
1045      INTEGER                         , INTENT(in   ) ::   nb , ndir
1046      !
1047      INTEGER  ::   ji, jj
1048      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
1049      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1050      !!---------------------------------------------------------------------- 
1051      !
1052      IF( before ) THEN
1053         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
1054      ELSE
1055         western_side  = (nb == 1).AND.(ndir == 1)
1056         eastern_side  = (nb == 1).AND.(ndir == 2)
1057         southern_side = (nb == 2).AND.(ndir == 1)
1058         northern_side = (nb == 2).AND.(ndir == 2)
1059         zrhoy = Agrif_Rhoy()
1060         zrhot = Agrif_rhot()
1061         ! Time indexes bounds for integration
1062         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1063         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1064         ! Polynomial interpolation coefficients:
1065         IF( bdy_tinterp == 1 ) THEN
1066            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1067               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1068         ELSEIF( bdy_tinterp == 2 ) THEN
1069            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1070               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
1071
1072         ELSE
1073            ztcoeff = 1
1074         ENDIF
1075         !   
1076         IF(western_side) THEN
1077            ubdy_w(j1:j2) = ubdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1078         ENDIF
1079         IF(eastern_side) THEN
1080            ubdy_e(j1:j2) = ubdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1081         ENDIF
1082         IF(southern_side) THEN
1083            ubdy_s(i1:i2) = ubdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1084         ENDIF
1085         IF(northern_side) THEN
1086            ubdy_n(i1:i2) = ubdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1087         ENDIF
1088         !           
1089         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1090            IF(western_side) THEN
1091               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
1092            ENDIF
1093            IF(eastern_side) THEN
1094               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
1095            ENDIF
1096            IF(southern_side) THEN
1097               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
1098            ENDIF
1099            IF(northern_side) THEN
1100               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
1101            ENDIF
1102         ENDIF
1103      ENDIF
1104      !
1105   END SUBROUTINE interpunb
1106
1107
1108   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
1109      !!----------------------------------------------------------------------
1110      !!                  ***  ROUTINE interpvnb  ***
1111      !!---------------------------------------------------------------------- 
1112      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1113      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1114      LOGICAL                         , INTENT(in   ) ::   before
1115      INTEGER                         , INTENT(in   ) ::   nb , ndir
1116      !
1117      INTEGER  ::   ji,jj
1118      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1119      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1120      !!---------------------------------------------------------------------- 
1121      !
1122      IF( before ) THEN
1123         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
1124      ELSE
1125         western_side  = (nb == 1).AND.(ndir == 1)
1126         eastern_side  = (nb == 1).AND.(ndir == 2)
1127         southern_side = (nb == 2).AND.(ndir == 1)
1128         northern_side = (nb == 2).AND.(ndir == 2)
1129         zrhox = Agrif_Rhox()
1130         zrhot = Agrif_rhot()
1131         ! Time indexes bounds for integration
1132         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1133         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1134         IF( bdy_tinterp == 1 ) THEN
1135            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1136               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1137         ELSEIF( bdy_tinterp == 2 ) THEN
1138            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1139               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
1140         ELSE
1141            ztcoeff = 1
1142         ENDIF
1143         !
1144         IF(western_side) THEN
1145            vbdy_w(j1:j2) = vbdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1146         ENDIF
1147         IF(eastern_side) THEN
1148            vbdy_e(j1:j2) = vbdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1149         ENDIF
1150         IF(southern_side) THEN
1151            vbdy_s(i1:i2) = vbdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1152         ENDIF
1153         IF(northern_side) THEN
1154            vbdy_n(i1:i2) = vbdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1155         ENDIF
1156         !           
1157         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1158            IF(western_side) THEN
1159               vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1160                     &                                  * vmask(i1,j1:j2,1)
1161            ENDIF
1162            IF(eastern_side) THEN
1163               vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1164                     &                                  * vmask(i1,j1:j2,1)
1165            ENDIF
1166            IF(southern_side) THEN
1167               vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1168                     &                                  * vmask(i1:i2,j1,1)
1169            ENDIF
1170            IF(northern_side) THEN
1171               vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1172                     &                                  * vmask(i1:i2,j1,1)
1173            ENDIF
1174         ENDIF
1175      ENDIF
1176      !
1177   END SUBROUTINE interpvnb
1178
1179
1180   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1181      !!----------------------------------------------------------------------
1182      !!                  ***  ROUTINE interpub2b  ***
1183      !!---------------------------------------------------------------------- 
1184      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1185      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1186      LOGICAL                         , INTENT(in   ) ::   before
1187      INTEGER                         , INTENT(in   ) ::   nb , ndir
1188      !
1189      INTEGER  ::   ji,jj
1190      REAL(wp) ::   zrhot, zt0, zt1,zat
1191      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1192      !!---------------------------------------------------------------------- 
1193      IF( before ) THEN
1194         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1195      ELSE
1196         western_side  = (nb == 1).AND.(ndir == 1)
1197         eastern_side  = (nb == 1).AND.(ndir == 2)
1198         southern_side = (nb == 2).AND.(ndir == 1)
1199         northern_side = (nb == 2).AND.(ndir == 2)
1200         zrhot = Agrif_rhot()
1201         ! Time indexes bounds for integration
1202         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1203         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1204         ! Polynomial interpolation coefficients:
1205         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1206            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1207         !
1208         IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1209         IF(eastern_side ) ubdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1210         IF(southern_side) ubdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1211         IF(northern_side) ubdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1212      ENDIF
1213      !
1214   END SUBROUTINE interpub2b
1215   
1216
1217   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1218      !!----------------------------------------------------------------------
1219      !!                  ***  ROUTINE interpvb2b  ***
1220      !!---------------------------------------------------------------------- 
1221      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1222      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1223      LOGICAL                         , INTENT(in   ) ::   before
1224      INTEGER                         , INTENT(in   ) ::   nb , ndir
1225      !
1226      INTEGER ::   ji,jj
1227      REAL(wp) ::   zrhot, zt0, zt1,zat
1228      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
1229      !!---------------------------------------------------------------------- 
1230      !
1231      IF( before ) THEN
1232         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1233      ELSE     
1234         western_side  = (nb == 1).AND.(ndir == 1)
1235         eastern_side  = (nb == 1).AND.(ndir == 2)
1236         southern_side = (nb == 2).AND.(ndir == 1)
1237         northern_side = (nb == 2).AND.(ndir == 2)
1238         zrhot = Agrif_rhot()
1239         ! Time indexes bounds for integration
1240         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1241         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1242         ! Polynomial interpolation coefficients:
1243         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1244            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1245         !
1246         IF(western_side )   vbdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1247         IF(eastern_side )   vbdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1248         IF(southern_side)   vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1249         IF(northern_side)   vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1250      ENDIF
1251      !     
1252   END SUBROUTINE interpvb2b
1253
1254
1255   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1256      !!----------------------------------------------------------------------
1257      !!                  ***  ROUTINE interpe3t  ***
1258      !!---------------------------------------------------------------------- 
1259      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1260      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1261      LOGICAL                              , INTENT(in   ) :: before
1262      INTEGER                              , INTENT(in   ) :: nb , ndir
1263      !
1264      INTEGER :: ji, jj, jk
1265      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1266      REAL(wp) :: ztmpmsk     
1267      !!---------------------------------------------------------------------- 
1268      !   
1269      IF( before ) THEN
1270         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1271      ELSE
1272         western_side  = (nb == 1).AND.(ndir == 1)
1273         eastern_side  = (nb == 1).AND.(ndir == 2)
1274         southern_side = (nb == 2).AND.(ndir == 1)
1275         northern_side = (nb == 2).AND.(ndir == 2)
1276
1277         DO jk = k1, k2
1278            DO jj = j1, j2
1279               DO ji = i1, i2
1280                  ! Get velocity mask at boundary edge points:
1281                  IF( western_side )   ztmpmsk = umask(ji    ,jj    ,1)
1282                  IF( eastern_side )   ztmpmsk = umask(nlci-2,jj    ,1)
1283                  IF( northern_side)   ztmpmsk = vmask(ji    ,nlcj-2,1)
1284                  IF( southern_side)   ztmpmsk = vmask(ji    ,2     ,1)
1285                  !
1286                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) )*ztmpmsk > 1.D-2) THEN
1287                     IF (western_side) THEN
1288                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1289                     ELSEIF (eastern_side) THEN
1290                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1291                     ELSEIF (southern_side) THEN
1292                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1293                     ELSEIF (northern_side) THEN
1294                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1295                     ENDIF
1296                     WRITE(numout,*) '      ptab(ji,jj,jk), e3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1297                     kindic_agr = kindic_agr + 1
1298                  ENDIF
1299               END DO
1300            END DO
1301         END DO
1302         !
1303      ENDIF
1304      !
1305   END SUBROUTINE interpe3t
1306
1307
1308   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1309      !!----------------------------------------------------------------------
1310      !!                  ***  ROUTINE interpumsk  ***
1311      !!---------------------------------------------------------------------- 
1312      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1313      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1314      LOGICAL                              , INTENT(in   ) ::   before
1315      INTEGER                              , INTENT(in   ) ::   nb , ndir
1316      !
1317      INTEGER ::   ji, jj, jk
1318      LOGICAL ::   western_side, eastern_side   
1319      !!---------------------------------------------------------------------- 
1320      !   
1321      IF( before ) THEN
1322         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1323      ELSE
1324         western_side = (nb == 1).AND.(ndir == 1)
1325         eastern_side = (nb == 1).AND.(ndir == 2)
1326         DO jk = k1, k2
1327            DO jj = j1, j2
1328               DO ji = i1, i2
1329                   ! Velocity mask at boundary edge points:
1330                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1331                     IF (western_side) THEN
1332                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1333                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1334                        kindic_agr = kindic_agr + 1
1335                     ELSEIF (eastern_side) THEN
1336                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1337                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1338                        kindic_agr = kindic_agr + 1
1339                     ENDIF
1340                  ENDIF
1341               END DO
1342            END DO
1343         END DO
1344         !
1345      ENDIF
1346      !
1347   END SUBROUTINE interpumsk
1348
1349
1350   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1351      !!----------------------------------------------------------------------
1352      !!                  ***  ROUTINE interpvmsk  ***
1353      !!---------------------------------------------------------------------- 
1354      INTEGER                              , INTENT(in   ) ::   i1,i2,j1,j2,k1,k2
1355      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1356      LOGICAL                              , INTENT(in   ) ::   before
1357      INTEGER                              , INTENT(in   ) :: nb , ndir
1358      !
1359      INTEGER ::   ji, jj, jk
1360      LOGICAL ::   northern_side, southern_side     
1361      !!---------------------------------------------------------------------- 
1362      !   
1363      IF( before ) THEN
1364         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1365      ELSE
1366         southern_side = (nb == 2).AND.(ndir == 1)
1367         northern_side = (nb == 2).AND.(ndir == 2)
1368         DO jk = k1, k2
1369            DO jj = j1, j2
1370               DO ji = i1, i2
1371                   ! Velocity mask at boundary edge points:
1372                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1373                     IF (southern_side) THEN
1374                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1375                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1376                        kindic_agr = kindic_agr + 1
1377                     ELSEIF (northern_side) THEN
1378                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1379                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1380                        kindic_agr = kindic_agr + 1
1381                     ENDIF
1382                  ENDIF
1383               END DO
1384            END DO
1385         END DO
1386         !
1387      ENDIF
1388      !
1389   END SUBROUTINE interpvmsk
1390
1391# if defined key_zdftke
1392
1393   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1394      !!----------------------------------------------------------------------
1395      !!                  ***  ROUTINE interavm  ***
1396      !!---------------------------------------------------------------------- 
1397      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
1398      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
1399      LOGICAL                              , INTENT(in   ) ::   before
1400      REAL(wp), DIMENSION(k1:k2) :: tabin
1401      REAL(wp) :: h_in(k1:k2)
1402      REAL(wp) :: h_out(1:jpk)
1403      REAL(wp) :: zrhoxy
1404      INTEGER  :: N_in, N_out, ji, jj, jk
1405      !!---------------------------------------------------------------------- 
1406      !     
1407      zrhoxy = Agrif_rhox()*Agrif_rhoy()
1408      IF (before) THEN         
1409         DO jk=k1,k2
1410            DO jj=j1,j2
1411              DO ji=i1,i2
1412                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1413              END DO
1414           END DO
1415        END DO
1416#ifdef key_vertical         
1417        DO jk=k1,k2
1418           DO jj=j1,j2
1419              DO ji=i1,i2
1420                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e1e2t(ji,jj) * e3w_n(ji,jj,jk) 
1421              END DO
1422           END DO
1423        END DO
1424#else
1425      ptab(i1:i2,j1:j2,k1:k2,2) = 0._wp
1426#endif
1427      ELSE 
1428#ifdef key_vertical         
1429         avm_k(i1:i2,j1:j2,1:jpk) = 0.
1430         DO jj=j1,j2
1431            DO ji=i1,i2
1432               N_in = 0
1433               DO jk=k1,k2 !k2 = jpk of parent grid
1434                  IF (ptab(ji,jj,jk,2) == 0) EXIT
1435                  N_in = N_in + 1
1436                  tabin(jk) = ptab(ji,jj,jk,1)
1437                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1e2t(ji,jj)*zrhoxy)
1438               END DO
1439               N_out = 0
1440               DO jk=1,jpk ! jpk of child grid
1441                  IF (wmask(ji,jj,jk) == 0) EXIT
1442                  N_out = N_out + 1
1443                  h_out(jk) = e3t_n(ji,jj,jk)
1444               ENDDO
1445               IF (N_in > 0) THEN
1446                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out)
1447               ENDIF
1448            ENDDO
1449         ENDDO
1450#else
1451         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1452#endif
1453      ENDIF
1454      !
1455   END SUBROUTINE interpavm
1456
1457# endif /* key_zdftke */
1458
1459#else
1460   !!----------------------------------------------------------------------
1461   !!   Empty module                                          no AGRIF zoom
1462   !!----------------------------------------------------------------------
1463CONTAINS
1464   SUBROUTINE Agrif_OPA_Interp_empty
1465      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
1466   END SUBROUTINE Agrif_OPA_Interp_empty
1467#endif
1468
1469   !!======================================================================
1470END MODULE agrif_opa_interp
Note: See TracBrowser for help on using the repository browser.