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 @ 8135

Last change on this file since 8135 was 8135, checked in by timgraham, 7 years ago

Added changes to code in NST_SRC from dev_r5803_UKMO_AGRIF_Vert_interp
  • Property svn:keywords set to Id
File size: 54.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      return
579     
580      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
581      IF( zalpha > 1. )   zalpha = 1.
582      !
583      Agrif_SpecialValue    = 0.e0
584      Agrif_UseSpecialValue = .TRUE.
585      !
586      CALL Agrif_Bc_variable(avm_id ,calledweight=zalpha, procname=interpavm)       
587      !
588      Agrif_UseSpecialValue = .FALSE.
589      !
590   END SUBROUTINE Agrif_tke
591   
592# endif
593
594   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
595      !!----------------------------------------------------------------------
596      !!   *** ROUTINE interptsn ***
597      !!----------------------------------------------------------------------
598      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
599      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
600      LOGICAL                                     , INTENT(in   ) ::   before
601      INTEGER                                     , INTENT(in   ) ::   nb , ndir
602      !
603      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
604      INTEGER  ::   imin, imax, jmin, jmax
605      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
606      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
607      LOGICAL :: western_side, eastern_side,northern_side,southern_side
608! VERTICAL REFINEMENT BEGIN
609      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
610      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
611      REAL(wp) :: h_in(k1:k2)
612      REAL(wp) :: h_out(1:jpk)
613      INTEGER :: N_in, N_out
614      REAL(wp) :: h_diff
615      REAL(wp) :: zrhoxy
616! VERTICAL REFINEMENT END
617
618      zrhoxy = Agrif_rhox()*Agrif_rhoy()
619      IF (before) THEN         
620         DO jn = n1,n2-1
621            DO jk=k1,k2
622               DO jj=j1,j2
623                  DO ji=i1,i2
624                     ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)
625                  END DO
626               END DO
627            END DO
628         END DO
629         DO jk=k1,k2
630            DO jj=j1,j2
631               DO ji=i1,i2
632                  ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
633               END DO
634            END DO
635         END DO
636
637      ELSE 
638! VERTICAL REFINEMENT BEGIN
639
640         ptab_child(:,:,:,:) = 0.
641         do jj=j1,j2
642         do ji=i1,i2
643           N_in = 0
644           DO jk=k1,k2 !k2 = jpk of parent grid
645             IF (ptab(ji,jj,jk,n2) == 0) EXIT
646             N_in = N_in + 1
647             tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2)
648             h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy)
649           END DO
650           N_out = 0
651           DO jk=1,jpk ! jpk of child grid
652             IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model
653             N_out = N_out + 1
654             h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
655           ENDDO
656           IF (N_in > 0) THEN
657             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
658!             if (h_diff > 0) then
659!               h_in(N_in+1) = h_diff
660!               N_in = N_in + 1
661!             else
662!               h_out(N_out+1) = -h_diff
663!               N_out = N_out + 1
664!             endif
665           ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) !what is this line for?????
666           do jn=1,jpts
667             call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
668           enddo
669!           if (abs(h_diff) > 1000.) then
670!           do jn=1,jpts
671!             do jk=1,N_out
672!             print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn)
673!             enddo
674!           enddo
675!         endif
676           ENDIF
677         enddo
678         enddo
679   
680
681! VERTICAL REFINEMENT END
682
683         !
684         western_side  = (nb == 1).AND.(ndir == 1)
685         eastern_side  = (nb == 1).AND.(ndir == 2)
686         southern_side = (nb == 2).AND.(ndir == 1)
687         northern_side = (nb == 2).AND.(ndir == 2)
688         !
689         zrhox = Agrif_Rhox()
690         !
691         zalpha1 = ( zrhox - 1. ) * 0.5
692         zalpha2 = 1. - zalpha1
693         !
694         zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. )
695         zalpha4 = 1. - zalpha3
696         !
697         zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
698         zalpha7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
699         zalpha5 = 1. - zalpha6 - zalpha7
700         !
701         imin = i1
702         imax = i2
703         jmin = j1
704         jmax = j2
705         !
706         ! Remove CORNERS
707         IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3
708         IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2
709         IF((nbondi == -1).OR.(nbondi == 2)) imin = 3
710         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2       
711         !
712         IF( eastern_side ) THEN
713            DO jn = 1, jpts
714               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)
715               DO jk = 1, jpkm1
716                  DO jj = jmin,jmax
717                     IF( umask(nlci-2,jj,jk) == 0._wp ) THEN
718                        tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
719                     ELSE
720                        tsa(nlci-1,jj,jk,jn)=(zalpha4*tsa(nlci,jj,jk,jn)+zalpha3*tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
721                        IF( un(nlci-2,jj,jk) > 0._wp ) THEN
722                           tsa(nlci-1,jj,jk,jn)=( zalpha6*tsa(nlci-2,jj,jk,jn)+zalpha5*tsa(nlci,jj,jk,jn) & 
723                                 + zalpha7*tsa(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk)
724                        ENDIF
725                     ENDIF
726                  END DO
727               END DO
728               tsa(nlci,j1:j2,k1:k2,jn) = 0._wp
729            END DO
730         ENDIF
731         !
732         IF( northern_side ) THEN           
733            DO jn = 1, jpts
734               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)
735               DO jk = 1, jpkm1
736                  DO ji = imin,imax
737                     IF( vmask(ji,nlcj-2,jk) == 0._wp ) THEN
738                        tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
739                     ELSE
740                        tsa(ji,nlcj-1,jk,jn)=(zalpha4*tsa(ji,nlcj,jk,jn)+zalpha3*tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)       
741                        IF (vn(ji,nlcj-2,jk) > 0._wp ) THEN
742                           tsa(ji,nlcj-1,jk,jn)=( zalpha6*tsa(ji,nlcj-2,jk,jn)+zalpha5*tsa(ji,nlcj,jk,jn)  &
743                                 + zalpha7*tsa(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk)
744                        ENDIF
745                     ENDIF
746                  END DO
747               END DO
748               tsa(i1:i2,nlcj,k1:k2,jn) = 0._wp
749            END DO
750         ENDIF
751         !
752         IF( western_side ) THEN           
753            DO jn = 1, jpts
754               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)
755               DO jk = 1, jpkm1
756                  DO jj = jmin,jmax
757                     IF( umask(2,jj,jk) == 0._wp ) THEN
758                        tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
759                     ELSE
760                        tsa(2,jj,jk,jn)=(zalpha4*tsa(1,jj,jk,jn)+zalpha3*tsa(3,jj,jk,jn))*tmask(2,jj,jk)       
761                        IF( un(2,jj,jk) < 0._wp ) THEN
762                           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)
763                        ENDIF
764                     ENDIF
765                  END DO
766               END DO
767               tsa(1,j1:j2,k1:k2,jn) = 0._wp
768            END DO
769         ENDIF
770         !
771         IF( southern_side ) THEN           
772            DO jn = 1, jpts
773               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)
774               DO jk = 1, jpk     
775                  DO ji=imin,imax
776                     IF( vmask(ji,2,jk) == 0._wp ) THEN
777                        tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
778                     ELSE
779                        tsa(ji,2,jk,jn)=(zalpha4*tsa(ji,1,jk,jn)+zalpha3*tsa(ji,3,jk,jn))*tmask(ji,2,jk)
780                        IF( vn(ji,2,jk) < 0._wp ) THEN
781                           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)
782                        ENDIF
783                     ENDIF
784                  END DO
785               END DO
786               tsa(i1:i2,1,k1:k2,jn) = 0._wp
787            END DO
788         ENDIF
789         !
790         ! Treatment of corners
791         !
792         ! East south
793         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
794            tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts)
795         ENDIF
796         ! East north
797         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
798            tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts)
799         ENDIF
800         ! West south
801         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
802            tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts)
803         ENDIF
804         ! West north
805         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
806            tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts)
807         ENDIF
808         !
809      ENDIF
810      !
811   END SUBROUTINE interptsn
812
813
814   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
815      !!----------------------------------------------------------------------
816      !!                  ***  ROUTINE interpsshn  ***
817      !!---------------------------------------------------------------------- 
818      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
819      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
820      LOGICAL                         , INTENT(in   ) ::   before
821      INTEGER                         , INTENT(in   ) ::   nb , ndir
822      !
823      LOGICAL :: western_side, eastern_side,northern_side,southern_side
824      !!---------------------------------------------------------------------- 
825      !
826      IF( before) THEN
827         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
828      ELSE
829         western_side  = (nb == 1).AND.(ndir == 1)
830         eastern_side  = (nb == 1).AND.(ndir == 2)
831         southern_side = (nb == 2).AND.(ndir == 1)
832         northern_side = (nb == 2).AND.(ndir == 2)
833         IF(western_side)  hbdy_w(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
834         IF(eastern_side)  hbdy_e(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
835         IF(southern_side) hbdy_s(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
836         IF(northern_side) hbdy_n(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
837      ENDIF
838      !
839   END SUBROUTINE interpsshn
840
841
842   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, before )
843      !!----------------------------------------------------------------------
844      !!   *** ROUTINE interpun ***
845      !!---------------------------------------------   
846      !!
847      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
848      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
849      LOGICAL, INTENT(in) :: before
850      INTEGER, INTENT(in) :: nb , ndir
851      !!
852      INTEGER :: ji,jj,jk
853      REAL(wp) :: zrhoy
854! VERTICAL REFINEMENT BEGIN
855      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
856      REAL(wp), DIMENSION(k1:k2) :: tabin
857      REAL(wp) :: h_in(k1:k2)
858      REAL(wp) :: h_out(1:jpk)
859      INTEGER :: N_in, N_out
860      REAL(wp) :: h_diff
861      LOGICAL :: western_side, eastern_side
862      INTEGER :: iref
863
864! VERTICAL REFINEMENT END
865      !!---------------------------------------------   
866      !
867      zrhoy = Agrif_rhoy()
868      IF (before) THEN 
869         !We can't use zero as the special value because we need to include zeros
870         !when interpolating the scale factors
871         IF(Agrif_UseSpecialValue) THEN
872             Agrif_SpecialValue = -999._wp
873         ELSE
874             Agrif_SpecialValue = 0._wp
875         ENDIF
876         DO jk=1,jpk
877            DO jj=j1,j2
878               DO ji=i1,i2
879                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) - &
880                                   & ((umask(ji,jj,jk)-1) * Agrif_SpecialValue)
881                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
882               END DO
883            END DO
884         END DO
885      ELSE
886! VERTICAL REFINEMENT BEGIN
887         western_side  = (nb == 1).AND.(ndir == 1)
888         eastern_side  = (nb == 1).AND.(ndir == 2)
889
890         Agrif_SpecialValue = 0._wp ! reset specialvalue to zero now interpolation completed
891
892         ptab_child(:,:,:) = 0.
893         DO jj=j1,j2
894            DO ji=i1,i2
895               iref = ji
896               IF (western_side) iref = 2
897               IF (eastern_side) iref = nlci-2
898
899               N_in = 0
900               DO jk=k1,k2
901                  IF (ptab(ji,jj,jk,2) == 0) EXIT
902                  N_in = N_in + 1
903                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
904                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj) * zrhoy) 
905              ENDDO
906         
907              IF (N_in == 0) THEN
908                 ptab_child(ji,jj,:) = 0.
909                 CYCLE
910              ENDIF
911         
912              N_out = 0
913              DO jk=1,jpk
914                 if (umask(ji,jj,jk) == 0) EXIT
915                 N_out = N_out + 1
916                 h_out(N_out) = e3u_n(ji,jj,jk)
917              ENDDO
918         
919              IF (N_out == 0) THEN
920                 ptab_child(ji,jj,:) = 0.
921                 CYCLE
922              ENDIF
923         
924!              IF (N_in * N_out > 0) THEN
925!                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
926! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
927!                 if (h_diff < 0.) then
928!                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
929!                    stop
930!                 endif
931!              ENDIF
932              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)
933         
934            ENDDO
935         ENDDO
936
937! in the following
938! remove division of ua by fs e3u (already done) and also zrhoy and e2u
939! VERTICAL REFINEMENT END
940
941         DO jk = 1, jpkm1
942            DO jj=j1,j2
943               ua(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
944!/(zrhoy*e2u(i1:i2,jj)))
945            END DO
946         END DO
947      ENDIF
948      !
949   END SUBROUTINE interpun
950
951
952   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, before )
953      !!----------------------------------------------------------------------
954      !!   *** ROUTINE interpvn ***
955      !!----------------------------------------------------------------------
956      !
957      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
958      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
959      LOGICAL, INTENT(in) :: before
960      INTEGER, INTENT(in) :: nb , ndir
961      !
962      INTEGER :: ji,jj,jk
963      REAL(wp) :: zrhox
964! VERTICAL REFINEMENT BEGIN
965      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
966      REAL(wp), DIMENSION(k1:k2) :: tabin
967      REAL(wp) :: h_in(k1:k2)
968      REAL(wp) :: h_out(1:jpk)
969      INTEGER :: N_in, N_out
970      REAL(wp) :: h_diff
971      LOGICAL :: northern_side,southern_side
972      INTEGER :: jref
973
974! VERTICAL REFINEMENT END
975      !!---------------------------------------------   
976      !     
977      zrhox = Agrif_rhox()
978      IF (before) THEN         
979         IF(Agrif_UseSpecialValue) THEN
980             Agrif_SpecialValue = -999._wp
981         ELSE
982             Agrif_SpecialValue = 0._wp
983         ENDIF
984         DO jk=k1,k2
985            DO jj=j1,j2
986               DO ji=i1,i2
987                  ptab(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
988                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
989               END DO
990            END DO
991         END DO
992      ELSE       
993! VERTICAL REFINEMENT BEGIN
994         ptab_child(:,:,:) = 0.
995         southern_side = (nb == 2).AND.(ndir == 1)
996         northern_side = (nb == 2).AND.(ndir == 2)
997
998         Agrif_SpecialValue = 0._wp !Reset special value to zero now interpolation is done
999
1000         do jj=j1,j2
1001            jref = jj
1002            IF (southern_side) jref = 2
1003            IF (northern_side) jref = nlcj-2
1004            do ji=i1,i2
1005
1006               N_in = 0
1007               do jk=k1,k2
1008                  if (ptab(ji,jj,jk,2) == 0) EXIT
1009                  N_in = N_in + 1
1010                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1011                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
1012               enddo
1013               IF (N_in == 0) THEN
1014                  ptab_child(ji,jj,:) = 0.
1015                  CYCLE
1016               ENDIF
1017         
1018               N_out = 0
1019               do jk=1,jpk
1020                  if (vmask(ji,jref,jk) == 0) EXIT
1021                  N_out = N_out + 1
1022                  h_out(N_out) = e3v_n(ji,jj,jk)
1023               enddo
1024               IF (N_out == 0) THEN
1025                 ptab_child(ji,jj,:) = 0.
1026                 CYCLE
1027               ENDIF
1028
1029!               IF (N_in * N_out > 0) THEN
1030!                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
1031! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
1032!                  if (h_diff < 0.) then
1033!                     print *,'CHECK YOUR BATHY interpvn...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
1034!                     stop
1035!                  endif
1036!               ENDIF
1037               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)
1038         
1039            enddo
1040         enddo
1041! in the following
1042! remove division of va by fs e3v, zrhox and e1v (already done)
1043! VERTICAL REFINEMENT END
1044         DO jk=1,jpkm1
1045            DO jj=j1,j2
1046               va(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
1047            END DO
1048         END DO
1049      ENDIF
1050      !       
1051   END SUBROUTINE interpvn
1052   
1053
1054   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
1055      !!----------------------------------------------------------------------
1056      !!                  ***  ROUTINE interpunb  ***
1057      !!---------------------------------------------------------------------- 
1058      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1059      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1060      LOGICAL                         , INTENT(in   ) ::   before
1061      INTEGER                         , INTENT(in   ) ::   nb , ndir
1062      !
1063      INTEGER  ::   ji, jj
1064      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
1065      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1066      !!---------------------------------------------------------------------- 
1067      !
1068      IF( before ) THEN
1069         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
1070      ELSE
1071         western_side  = (nb == 1).AND.(ndir == 1)
1072         eastern_side  = (nb == 1).AND.(ndir == 2)
1073         southern_side = (nb == 2).AND.(ndir == 1)
1074         northern_side = (nb == 2).AND.(ndir == 2)
1075         zrhoy = Agrif_Rhoy()
1076         zrhot = Agrif_rhot()
1077         ! Time indexes bounds for integration
1078         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1079         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1080         ! Polynomial interpolation coefficients:
1081         IF( bdy_tinterp == 1 ) THEN
1082            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1083               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1084         ELSEIF( bdy_tinterp == 2 ) THEN
1085            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1086               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
1087
1088         ELSE
1089            ztcoeff = 1
1090         ENDIF
1091         !   
1092         IF(western_side) THEN
1093            ubdy_w(j1:j2) = ubdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1094         ENDIF
1095         IF(eastern_side) THEN
1096            ubdy_e(j1:j2) = ubdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1097         ENDIF
1098         IF(southern_side) THEN
1099            ubdy_s(i1:i2) = ubdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1100         ENDIF
1101         IF(northern_side) THEN
1102            ubdy_n(i1:i2) = ubdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1103         ENDIF
1104         !           
1105         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1106            IF(western_side) THEN
1107               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
1108            ENDIF
1109            IF(eastern_side) THEN
1110               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
1111            ENDIF
1112            IF(southern_side) THEN
1113               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
1114            ENDIF
1115            IF(northern_side) THEN
1116               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
1117            ENDIF
1118         ENDIF
1119      ENDIF
1120      !
1121   END SUBROUTINE interpunb
1122
1123
1124   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
1125      !!----------------------------------------------------------------------
1126      !!                  ***  ROUTINE interpvnb  ***
1127      !!---------------------------------------------------------------------- 
1128      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1129      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1130      LOGICAL                         , INTENT(in   ) ::   before
1131      INTEGER                         , INTENT(in   ) ::   nb , ndir
1132      !
1133      INTEGER  ::   ji,jj
1134      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
1135      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1136      !!---------------------------------------------------------------------- 
1137      !
1138      IF( before ) THEN
1139         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
1140      ELSE
1141         western_side  = (nb == 1).AND.(ndir == 1)
1142         eastern_side  = (nb == 1).AND.(ndir == 2)
1143         southern_side = (nb == 2).AND.(ndir == 1)
1144         northern_side = (nb == 2).AND.(ndir == 2)
1145         zrhox = Agrif_Rhox()
1146         zrhot = Agrif_rhot()
1147         ! Time indexes bounds for integration
1148         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1149         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
1150         IF( bdy_tinterp == 1 ) THEN
1151            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1152               &               - zt0**2._wp * (       zt0 - 1._wp)        )
1153         ELSEIF( bdy_tinterp == 2 ) THEN
1154            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1155               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
1156         ELSE
1157            ztcoeff = 1
1158         ENDIF
1159         !
1160         IF(western_side) THEN
1161            vbdy_w(j1:j2) = vbdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1162         ENDIF
1163         IF(eastern_side) THEN
1164            vbdy_e(j1:j2) = vbdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2) 
1165         ENDIF
1166         IF(southern_side) THEN
1167            vbdy_s(i1:i2) = vbdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1168         ENDIF
1169         IF(northern_side) THEN
1170            vbdy_n(i1:i2) = vbdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1) 
1171         ENDIF
1172         !           
1173         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1174            IF(western_side) THEN
1175               vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1176                     &                                  * vmask(i1,j1:j2,1)
1177            ENDIF
1178            IF(eastern_side) THEN
1179               vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
1180                     &                                  * vmask(i1,j1:j2,1)
1181            ENDIF
1182            IF(southern_side) THEN
1183               vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1184                     &                                  * vmask(i1:i2,j1,1)
1185            ENDIF
1186            IF(northern_side) THEN
1187               vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
1188                     &                                  * vmask(i1:i2,j1,1)
1189            ENDIF
1190         ENDIF
1191      ENDIF
1192      !
1193   END SUBROUTINE interpvnb
1194
1195
1196   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1197      !!----------------------------------------------------------------------
1198      !!                  ***  ROUTINE interpub2b  ***
1199      !!---------------------------------------------------------------------- 
1200      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1201      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1202      LOGICAL                         , INTENT(in   ) ::   before
1203      INTEGER                         , INTENT(in   ) ::   nb , ndir
1204      !
1205      INTEGER  ::   ji,jj
1206      REAL(wp) ::   zrhot, zt0, zt1,zat
1207      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
1208      !!---------------------------------------------------------------------- 
1209      IF( before ) THEN
1210         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1211      ELSE
1212         western_side  = (nb == 1).AND.(ndir == 1)
1213         eastern_side  = (nb == 1).AND.(ndir == 2)
1214         southern_side = (nb == 2).AND.(ndir == 1)
1215         northern_side = (nb == 2).AND.(ndir == 2)
1216         zrhot = Agrif_rhot()
1217         ! Time indexes bounds for integration
1218         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1219         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1220         ! Polynomial interpolation coefficients:
1221         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1222            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1223         !
1224         IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1225         IF(eastern_side ) ubdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1226         IF(southern_side) ubdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1227         IF(northern_side) ubdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1228      ENDIF
1229      !
1230   END SUBROUTINE interpub2b
1231   
1232
1233   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1234      !!----------------------------------------------------------------------
1235      !!                  ***  ROUTINE interpvb2b  ***
1236      !!---------------------------------------------------------------------- 
1237      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1238      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1239      LOGICAL                         , INTENT(in   ) ::   before
1240      INTEGER                         , INTENT(in   ) ::   nb , ndir
1241      !
1242      INTEGER ::   ji,jj
1243      REAL(wp) ::   zrhot, zt0, zt1,zat
1244      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
1245      !!---------------------------------------------------------------------- 
1246      !
1247      IF( before ) THEN
1248         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1249      ELSE     
1250         western_side  = (nb == 1).AND.(ndir == 1)
1251         eastern_side  = (nb == 1).AND.(ndir == 2)
1252         southern_side = (nb == 2).AND.(ndir == 1)
1253         northern_side = (nb == 2).AND.(ndir == 2)
1254         zrhot = Agrif_rhot()
1255         ! Time indexes bounds for integration
1256         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1257         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1258         ! Polynomial interpolation coefficients:
1259         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1260            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
1261         !
1262         IF(western_side )   vbdy_w(j1:j2) = zat * ptab(i1,j1:j2) 
1263         IF(eastern_side )   vbdy_e(j1:j2) = zat * ptab(i1,j1:j2) 
1264         IF(southern_side)   vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
1265         IF(northern_side)   vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
1266      ENDIF
1267      !     
1268   END SUBROUTINE interpvb2b
1269
1270
1271   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1272      !!----------------------------------------------------------------------
1273      !!                  ***  ROUTINE interpe3t  ***
1274      !!---------------------------------------------------------------------- 
1275      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
1276      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1277      LOGICAL                              , INTENT(in   ) :: before
1278      INTEGER                              , INTENT(in   ) :: nb , ndir
1279      !
1280      INTEGER :: ji, jj, jk
1281      LOGICAL :: western_side, eastern_side, northern_side, southern_side
1282      REAL(wp) :: ztmpmsk     
1283      !!---------------------------------------------------------------------- 
1284      !   
1285      IF( before ) THEN
1286         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1287      ELSE
1288         western_side  = (nb == 1).AND.(ndir == 1)
1289         eastern_side  = (nb == 1).AND.(ndir == 2)
1290         southern_side = (nb == 2).AND.(ndir == 1)
1291         northern_side = (nb == 2).AND.(ndir == 2)
1292
1293         DO jk = k1, k2
1294            DO jj = j1, j2
1295               DO ji = i1, i2
1296                  ! Get velocity mask at boundary edge points:
1297                  IF( western_side )   ztmpmsk = umask(ji    ,jj    ,1)
1298                  IF( eastern_side )   ztmpmsk = umask(nlci-2,jj    ,1)
1299                  IF( northern_side)   ztmpmsk = vmask(ji    ,nlcj-2,1)
1300                  IF( southern_side)   ztmpmsk = vmask(ji    ,2     ,1)
1301                  !
1302                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) )*ztmpmsk > 1.D-2) THEN
1303                     IF (western_side) THEN
1304                        WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1305                     ELSEIF (eastern_side) THEN
1306                        WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1307                     ELSEIF (southern_side) THEN
1308                        WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1309                     ELSEIF (northern_side) THEN
1310                        WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1311                     ENDIF
1312                     WRITE(numout,*) '      ptab(ji,jj,jk), e3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1313                     kindic_agr = kindic_agr + 1
1314                  ENDIF
1315               END DO
1316            END DO
1317         END DO
1318         !
1319      ENDIF
1320      !
1321   END SUBROUTINE interpe3t
1322
1323
1324   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1325      !!----------------------------------------------------------------------
1326      !!                  ***  ROUTINE interpumsk  ***
1327      !!---------------------------------------------------------------------- 
1328      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1329      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1330      LOGICAL                              , INTENT(in   ) ::   before
1331      INTEGER                              , INTENT(in   ) ::   nb , ndir
1332      !
1333      INTEGER ::   ji, jj, jk
1334      LOGICAL ::   western_side, eastern_side   
1335      !!---------------------------------------------------------------------- 
1336      !   
1337      IF( before ) THEN
1338         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1339      ELSE
1340         western_side = (nb == 1).AND.(ndir == 1)
1341         eastern_side = (nb == 1).AND.(ndir == 2)
1342         DO jk = k1, k2
1343            DO jj = j1, j2
1344               DO ji = i1, i2
1345                   ! Velocity mask at boundary edge points:
1346                  IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1347                     IF (western_side) THEN
1348                        WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1349                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1350                        kindic_agr = kindic_agr + 1
1351                     ELSEIF (eastern_side) THEN
1352                        WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1353                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1354                        kindic_agr = kindic_agr + 1
1355                     ENDIF
1356                  ENDIF
1357               END DO
1358            END DO
1359         END DO
1360         !
1361      ENDIF
1362      !
1363   END SUBROUTINE interpumsk
1364
1365
1366   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1367      !!----------------------------------------------------------------------
1368      !!                  ***  ROUTINE interpvmsk  ***
1369      !!---------------------------------------------------------------------- 
1370      INTEGER                              , INTENT(in   ) ::   i1,i2,j1,j2,k1,k2
1371      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1372      LOGICAL                              , INTENT(in   ) ::   before
1373      INTEGER                              , INTENT(in   ) :: nb , ndir
1374      !
1375      INTEGER ::   ji, jj, jk
1376      LOGICAL ::   northern_side, southern_side     
1377      !!---------------------------------------------------------------------- 
1378      !   
1379      IF( before ) THEN
1380         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1381      ELSE
1382         southern_side = (nb == 2).AND.(ndir == 1)
1383         northern_side = (nb == 2).AND.(ndir == 2)
1384         DO jk = k1, k2
1385            DO jj = j1, j2
1386               DO ji = i1, i2
1387                   ! Velocity mask at boundary edge points:
1388                  IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1389                     IF (southern_side) THEN
1390                        WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1391                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1392                        kindic_agr = kindic_agr + 1
1393                     ELSEIF (northern_side) THEN
1394                        WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1395                        WRITE(numout,*) '      masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1396                        kindic_agr = kindic_agr + 1
1397                     ENDIF
1398                  ENDIF
1399               END DO
1400            END DO
1401         END DO
1402         !
1403      ENDIF
1404      !
1405   END SUBROUTINE interpvmsk
1406
1407# if defined key_zdftke
1408
1409   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )
1410      !!----------------------------------------------------------------------
1411      !!                  ***  ROUTINE interavm  ***
1412      !!---------------------------------------------------------------------- 
1413      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1414      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1415      LOGICAL                              , INTENT(in   ) ::   before
1416      !!---------------------------------------------------------------------- 
1417      !     
1418      IF( before ) THEN
1419         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1420      ELSE
1421         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1422      ENDIF
1423      !
1424   END SUBROUTINE interpavm
1425
1426# endif /* key_zdftke */
1427
1428#else
1429   !!----------------------------------------------------------------------
1430   !!   Empty module                                          no AGRIF zoom
1431   !!----------------------------------------------------------------------
1432CONTAINS
1433   SUBROUTINE Agrif_OPA_Interp_empty
1434      WRITE(*,*)  'agrif_opa_interp : You should not have seen this print! error?'
1435   END SUBROUTINE Agrif_OPA_Interp_empty
1436#endif
1437
1438   !!======================================================================
1439END MODULE agrif_opa_interp
Note: See TracBrowser for help on using the repository browser.