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_oce_interp.F90 in NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90 @ 12152

Last change on this file since 12152 was 12152, checked in by jchanut, 4 years ago

#2222: fixes linear vertical interpolation of viscosities

  • Property svn:keywords set to Id
File size: 44.8 KB
RevLine 
[9570]1MODULE agrif_oce_interp
[1605]2   !!======================================================================
[9570]3   !!                   ***  MODULE  agrif_oce_interp  ***
[9019]4   !! AGRIF: interpolation package for the ocean dynamics (OPA)
[1605]5   !!======================================================================
[9019]6   !! History :  2.0  !  2002-06  (L. Debreu)  Original cade
[1605]7   !!            3.2  !  2009-04  (R. Benshila)
[5656]8   !!            3.6  !  2014-09  (R. Benshila)
[1605]9   !!----------------------------------------------------------------------
[7646]10#if defined key_agrif
[1605]11   !!----------------------------------------------------------------------
12   !!   'key_agrif'                                              AGRIF zoom
13   !!----------------------------------------------------------------------
14   !!   Agrif_tra     :
15   !!   Agrif_dyn     :
[9019]16   !!   Agrif_ssh     :
17   !!   Agrif_dyn_ts  :
18   !!   Agrif_dta_ts  :
19   !!   Agrif_ssh_ts  :
20   !!   Agrif_avm     :
[1605]21   !!   interpu       :
22   !!   interpv       :
23   !!----------------------------------------------------------------------
[636]24   USE par_oce
25   USE oce
26   USE dom_oce     
[6140]27   USE zdf_oce
[782]28   USE agrif_oce
[1605]29   USE phycst
[9031]30   USE dynspg_ts, ONLY: un_adv, vn_adv
[6140]31   !
[1605]32   USE in_out_manager
[9570]33   USE agrif_oce_sponge
[2715]34   USE lib_mpp
[11590]35   USE vremap
[5656]36 
[636]37   IMPLICIT NONE
38   PRIVATE
[4292]39
[11574]40   PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_dyn_ts_flux, Agrif_ssh_ts, Agrif_dta_ts
[9057]41   PUBLIC   Agrif_tra, Agrif_avm
[9019]42   PUBLIC   interpun , interpvn
[9057]43   PUBLIC   interptsn, interpsshn, interpavm
[9019]44   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
[12119]45   PUBLIC   interpe3t
[11741]46#if defined key_vertical
47   PUBLIC   interpht0, interpmbkt
48# endif
[11603]49   INTEGER ::   bdy_tinterp = 0
50
[1605]51#  include "vectopt_loop_substitute.h90"
[1156]52   !!----------------------------------------------------------------------
[9598]53   !! NEMO/NST 4.0 , NEMO Consortium (2018)
[1156]54   !! $Id$
[10068]55   !! Software governed by the CeCILL license (see ./LICENSE)
[1156]56   !!----------------------------------------------------------------------
[5656]57CONTAINS
58
[782]59   SUBROUTINE Agrif_tra
[1605]60      !!----------------------------------------------------------------------
[5656]61      !!                  ***  ROUTINE Agrif_tra  ***
[1605]62      !!----------------------------------------------------------------------
[636]63      !
[1605]64      IF( Agrif_Root() )   RETURN
[6140]65      !
66      Agrif_SpecialValue    = 0._wp
[636]67      Agrif_UseSpecialValue = .TRUE.
[6140]68      !
[5656]69      CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
[6140]70      !
[636]71      Agrif_UseSpecialValue = .FALSE.
[1605]72      !
[636]73   END SUBROUTINE Agrif_tra
74
[1605]75
[636]76   SUBROUTINE Agrif_dyn( kt )
[1605]77      !!----------------------------------------------------------------------
78      !!                  ***  ROUTINE Agrif_DYN  ***
79      !!---------------------------------------------------------------------- 
80      INTEGER, INTENT(in) ::   kt
[6140]81      !
82      INTEGER ::   ji, jj, jk       ! dummy loop indices
[9057]83      INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2
[9019]84      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb
[1605]85      !!---------------------------------------------------------------------- 
[6140]86      !
[1605]87      IF( Agrif_Root() )   RETURN
[6140]88      !
89      Agrif_SpecialValue    = 0._wp
[5656]90      Agrif_UseSpecialValue = ln_spc_dyn
[6140]91      !
92      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
93      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
94      !
[5656]95      Agrif_UseSpecialValue = .FALSE.
[6140]96      !
[11574]97      ! --- West --- !
98      ibdy1 = 2
99      ibdy2 = 1+nbghostcells 
100      !
101      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
102         DO ji = mi0(ibdy1), mi1(ibdy2)
103            ua_b(ji,:) = 0._wp
[782]104
[6140]105            DO jk = 1, jpkm1
106               DO jj = 1, jpj
[11574]107                  ua_b(ji,jj) = ua_b(ji,jj) + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
[5930]108               END DO
[636]109            END DO
[11574]110
[6140]111            DO jj = 1, jpj
[11574]112               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj)
[636]113            END DO
[11574]114         END DO
115      ENDIF
116      !
117      DO ji = mi0(ibdy1), mi1(ibdy2)
118         zub(ji,:) = 0._wp    ! Correct transport
[9057]119         DO jk = 1, jpkm1
120            DO jj = 1, jpj
[11574]121               zub(ji,jj) = zub(ji,jj) & 
122                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk)*umask(ji,jj,jk)
[9057]123            END DO
124         END DO
125         DO jj=1,jpj
[11574]126            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
[9057]127         END DO
128           
129         DO jk = 1, jpkm1
130            DO jj = 1, jpj
[11574]131               ua(ji,jj,jk) = ( ua(ji,jj,jk) + ua_b(ji,jj)-zub(ji,jj)) * umask(ji,jj,jk)
[9057]132            END DO
133         END DO
[11574]134      END DO
[9057]135           
[11574]136      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
137         DO ji = mi0(ibdy1), mi1(ibdy2)
138            zvb(ji,:) = 0._wp
[9019]139            DO jk = 1, jpkm1
140               DO jj = 1, jpj
[11574]141                  zvb(ji,jj) = zvb(ji,jj) + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
[9019]142               END DO
[636]143            END DO
[9057]144            DO jj = 1, jpj
[11574]145               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
[636]146            END DO
[6140]147            DO jk = 1, jpkm1
148               DO jj = 1, jpj
[11574]149                  va(ji,jj,jk) = ( va(ji,jj,jk) + va_b(ji,jj)-zvb(ji,jj))*vmask(ji,jj,jk)
[5930]150               END DO
151            END DO
[11574]152         END DO
[636]153      ENDIF
[390]154
[9019]155      ! --- East --- !
[11574]156      ibdy1 = jpiglo-1-nbghostcells
157      ibdy2 = jpiglo-2 
158      !
159      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
160         DO ji = mi0(ibdy1), mi1(ibdy2)
161            ua_b(ji,:) = 0._wp
[9057]162            DO jk = 1, jpkm1
163               DO jj = 1, jpj
[11574]164                  ua_b(ji,jj) = ua_b(ji,jj) & 
165                      & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
[5930]166               END DO
[636]167            END DO
[9057]168            DO jj = 1, jpj
[11574]169               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj)
[5930]170            END DO
[11574]171         END DO
172      ENDIF
173      !
174      DO ji = mi0(ibdy1), mi1(ibdy2)
175         zub(ji,:) = 0._wp    ! Correct transport
[9031]176         DO jk = 1, jpkm1
177            DO jj = 1, jpj
[11574]178               zub(ji,jj) = zub(ji,jj) & 
179                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk) * umask(ji,jj,jk)
[9031]180            END DO
181         END DO
[9057]182         DO jj=1,jpj
[11574]183            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
[9031]184         END DO
[9057]185           
[9031]186         DO jk = 1, jpkm1
187            DO jj = 1, jpj
[11574]188               ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
189                 & + ua_b(ji,jj)-zub(ji,jj))*umask(ji,jj,jk)
[9031]190            END DO
191         END DO
[11574]192      END DO
[9057]193           
[11574]194      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
195         ibdy1 = jpiglo-nbghostcells
196         ibdy2 = jpiglo-1 
197         DO ji = mi0(ibdy1), mi1(ibdy2)
198            zvb(ji,:) = 0._wp
[6140]199            DO jk = 1, jpkm1
200               DO jj = 1, jpj
[11574]201                  zvb(ji,jj) = zvb(ji,jj) &
202                     & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
[5930]203               END DO
204            END DO
[9019]205            DO jj = 1, jpj
[11574]206               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
[4486]207            END DO
[6140]208            DO jk = 1, jpkm1
209               DO jj = 1, jpj
[11574]210                  va(ji,jj,jk) = ( va(ji,jj,jk) & 
211                      & + va_b(ji,jj)-zvb(ji,jj)) * vmask(ji,jj,jk)
[5930]212               END DO
213            END DO
[11574]214         END DO
[636]215      ENDIF
[390]216
[9019]217      ! --- South --- !
[11574]218      jbdy1 = 2
219      jbdy2 = 1+nbghostcells 
220      !
221      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
222         DO jj = mj0(jbdy1), mj1(jbdy2)
223            va_b(:,jj) = 0._wp
[6140]224            DO jk = 1, jpkm1
225               DO ji = 1, jpi
[11574]226                  va_b(ji,jj) = va_b(ji,jj) & 
227                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
[5930]228               END DO
[636]229            END DO
230            DO ji=1,jpi
[11574]231               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)     
[636]232            END DO
[11574]233         END DO
234      ENDIF
235      !
236      DO jj = mj0(jbdy1), mj1(jbdy2)
237         zvb(:,jj) = 0._wp    ! Correct transport
[9031]238         DO jk=1,jpkm1
239            DO ji=1,jpi
[11574]240               zvb(ji,jj) = zvb(ji,jj) & 
241                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
[636]242            END DO
[9057]243         END DO
244         DO ji = 1, jpi
[11574]245            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
[9057]246         END DO
[9134]247
[9057]248         DO jk = 1, jpkm1
[6140]249            DO ji = 1, jpi
[11574]250               va(ji,jj,jk) = ( va(ji,jj,jk) & 
251                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk)
[636]252            END DO
[9057]253         END DO
[11574]254      END DO
[9057]255           
[11574]256      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
257         DO jj = mj0(jbdy1), mj1(jbdy2)
258            zub(:,jj) = 0._wp
[6140]259            DO jk = 1, jpkm1
260               DO ji = 1, jpi
[11574]261                  zub(ji,jj) = zub(ji,jj) & 
262                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
[5930]263               END DO
264            END DO
[9057]265            DO ji = 1, jpi
[11574]266               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
[9057]267            END DO
268               
269            DO jk = 1, jpkm1
[6140]270               DO ji = 1, jpi
[11574]271                  ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
272                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk)
[5930]273               END DO
[9057]274            END DO
[11574]275         END DO
[636]276      ENDIF
[390]277
[9019]278      ! --- North --- !
[11574]279      jbdy1 = jpjglo-1-nbghostcells
280      jbdy2 = jpjglo-2 
281      !
282      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
283         DO jj = mj0(jbdy1), mj1(jbdy2)
284            va_b(:,jj) = 0._wp
[6140]285            DO jk = 1, jpkm1
286               DO ji = 1, jpi
[11574]287                  va_b(ji,jj) = va_b(ji,jj) & 
288                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
[5930]289               END DO
[636]290            END DO
[9057]291            DO ji=1,jpi
[11574]292               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)
[636]293            END DO
[11574]294         END DO
295      ENDIF
296      !
297      DO jj = mj0(jbdy1), mj1(jbdy2)
298         zvb(:,jj) = 0._wp    ! Correct transport
[9057]299         DO jk=1,jpkm1
300            DO ji=1,jpi
[11574]301               zvb(ji,jj) = zvb(ji,jj) & 
302                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
[9057]303            END DO
304         END DO
305         DO ji = 1, jpi
[11574]306            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
[9057]307         END DO
[9134]308
[9031]309         DO jk = 1, jpkm1
310            DO ji = 1, jpi
[11574]311               va(ji,jj,jk) = ( va(ji,jj,jk) & 
312                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk)
[636]313            END DO
[9057]314         END DO
[11574]315      END DO
[9057]316           
[11574]317      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
318         jbdy1 = jpjglo-nbghostcells
319         jbdy2 = jpjglo-1
320         DO jj = mj0(jbdy1), mj1(jbdy2)
321            zub(:,jj) = 0._wp
[9057]322            DO jk = 1, jpkm1
323               DO ji = 1, jpi
[11574]324                  zub(ji,jj) = zub(ji,jj) & 
325                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
[9057]326               END DO
327            END DO
[6140]328            DO ji = 1, jpi
[11574]329               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
[636]330            END DO
[9057]331               
[6140]332            DO jk = 1, jpkm1
333               DO ji = 1, jpi
[11574]334                  ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
335                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk)
[5930]336               END DO
337            END DO
[11574]338         END DO
[636]339      ENDIF
[2715]340      !
[636]341   END SUBROUTINE Agrif_dyn
[390]342
[6140]343
[4486]344   SUBROUTINE Agrif_dyn_ts( jn )
[4292]345      !!----------------------------------------------------------------------
346      !!                  ***  ROUTINE Agrif_dyn_ts  ***
347      !!---------------------------------------------------------------------- 
[4486]348      INTEGER, INTENT(in) ::   jn
[4292]349      !!
350      INTEGER :: ji, jj
[11574]351      INTEGER :: istart, iend, jstart, jend
[4486]352      !!---------------------------------------------------------------------- 
[6140]353      !
[4486]354      IF( Agrif_Root() )   RETURN
[9057]355      !
[11574]356      !--- West ---!
357      istart = 2
358      iend   = nbghostcells+1
359      DO ji = mi0(istart), mi1(iend)
[4486]360         DO jj=1,jpj
[11574]361            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
362            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
[4486]363         END DO
[11574]364      END DO
[6140]365      !
[11574]366      !--- East ---!
367      istart = jpiglo-nbghostcells
368      iend   = jpiglo-1
369      DO ji = mi0(istart), mi1(iend)
[4486]370         DO jj=1,jpj
[11574]371            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
[4486]372         END DO
[11574]373      END DO
374      istart = jpiglo-nbghostcells-1
375      iend   = jpiglo-2
376      DO ji = mi0(istart), mi1(iend)
377         DO jj=1,jpj
378            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
379         END DO
380      END DO
[6140]381      !
[11574]382      !--- South ---!
383      jstart = 2
384      jend   = nbghostcells+1
385      DO jj = mj0(jstart), mj1(jend)
[4486]386         DO ji=1,jpi
[11574]387            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
388            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
[4486]389         END DO
[11574]390      END DO
[6140]391      !
[11574]392      !--- North ---!
393      jstart = jpjglo-nbghostcells
394      jend   = jpjglo-1
395      DO jj = mj0(jstart), mj1(jend)
[4486]396         DO ji=1,jpi
[11574]397            ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
[4486]398         END DO
[11574]399      END DO
400      jstart = jpjglo-nbghostcells-1
401      jend   = jpjglo-2
402      DO jj = mj0(jstart), mj1(jend)
403         DO ji=1,jpi
404            va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
405         END DO
406      END DO
[4486]407      !
408   END SUBROUTINE Agrif_dyn_ts
409
[11574]410   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv )
411      !!----------------------------------------------------------------------
412      !!                  ***  ROUTINE Agrif_dyn_ts_flux  ***
413      !!---------------------------------------------------------------------- 
414      INTEGER, INTENT(in) ::   jn
415      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv
416      !!
417      INTEGER :: ji, jj
418      INTEGER :: istart, iend, jstart, jend
419      !!---------------------------------------------------------------------- 
420      !
421      IF( Agrif_Root() )   RETURN
422      !
423      !--- West ---!
424      istart = 2
425      iend   = nbghostcells+1
426      DO ji = mi0(istart), mi1(iend)
427         DO jj=1,jpj
428            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
429            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
430         END DO
431      END DO
432      !
433      !--- East ---!
434      istart = jpiglo-nbghostcells
435      iend   = jpiglo-1
436      DO ji = mi0(istart), mi1(iend)
437         DO jj=1,jpj
438            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
439         END DO
440      END DO
441      istart = jpiglo-nbghostcells-1
442      iend   = jpiglo-2
443      DO ji = mi0(istart), mi1(iend)
444         DO jj=1,jpj
445            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
446         END DO
447      END DO
448      !
449      !--- South ---!
450      jstart = 2
451      jend   = nbghostcells+1
452      DO jj = mj0(jstart), mj1(jend)
453         DO ji=1,jpi
454            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
455            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
456         END DO
457      END DO
458      !
459      !--- North ---!
460      jstart = jpjglo-nbghostcells
461      jend   = jpjglo-1
462      DO jj = mj0(jstart), mj1(jend)
463         DO ji=1,jpi
464            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
465         END DO
466      END DO
467      jstart = jpjglo-nbghostcells-1
468      jend   = jpjglo-2
469      DO jj = mj0(jstart), mj1(jend)
470         DO ji=1,jpi
471            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
472         END DO
473      END DO
474      !
475   END SUBROUTINE Agrif_dyn_ts_flux
[6140]476
[4486]477   SUBROUTINE Agrif_dta_ts( kt )
478      !!----------------------------------------------------------------------
479      !!                  ***  ROUTINE Agrif_dta_ts  ***
480      !!---------------------------------------------------------------------- 
481      INTEGER, INTENT(in) ::   kt
482      !!
483      INTEGER :: ji, jj
484      LOGICAL :: ll_int_cons
[4292]485      !!---------------------------------------------------------------------- 
[6140]486      !
[4292]487      IF( Agrif_Root() )   RETURN
[6140]488      !
489      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
490      !
[9031]491      ! Enforce volume conservation if no time refinement: 
492      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE. 
[6140]493      !
[4486]494      ! Interpolate barotropic fluxes
[11574]495      Agrif_SpecialValue = 0._wp
[4486]496      Agrif_UseSpecialValue = ln_spc_dyn
[6140]497      !
[11574]498      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
499      utint_stage(:,:) = 0
500      vtint_stage(:,:) = 0
501      !
[6140]502      IF( ll_int_cons ) THEN  ! Conservative interpolation
[9019]503         ! order matters here !!!!!!
[6140]504         CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
[11574]505         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 
506         !
[11603]507         bdy_tinterp = 1
[6140]508         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
[11603]509         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  ) 
[11574]510         !
[11603]511         bdy_tinterp = 2
[6140]512         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
[11603]513         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )   
[4486]514      ELSE ! Linear interpolation
[11574]515         !
516         ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp 
[9031]517         CALL Agrif_Bc_variable( unb_id, procname=interpunb )
518         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb )
[4486]519      ENDIF
520      Agrif_UseSpecialValue = .FALSE.
[5656]521      !
[4486]522   END SUBROUTINE Agrif_dta_ts
523
[6140]524
[2486]525   SUBROUTINE Agrif_ssh( kt )
526      !!----------------------------------------------------------------------
[9031]527      !!                  ***  ROUTINE Agrif_ssh  ***
[2486]528      !!---------------------------------------------------------------------- 
529      INTEGER, INTENT(in) ::   kt
[9019]530      !
[11574]531      INTEGER  :: ji, jj
532      INTEGER  :: istart, iend, jstart, jend
[2486]533      !!---------------------------------------------------------------------- 
[6140]534      !
[2486]535      IF( Agrif_Root() )   RETURN
[9031]536      !     
[9116]537      ! Linear time interpolation of sea level
[9031]538      !
539      Agrif_SpecialValue    = 0._wp
540      Agrif_UseSpecialValue = .TRUE.
541      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
542      Agrif_UseSpecialValue = .FALSE.
543      !
[9116]544      ! --- West --- !
[11574]545      istart = 2
546      iend   = 1 + nbghostcells
547      DO ji = mi0(istart), mi1(iend)
[9019]548         DO jj = 1, jpj
[11574]549            ssha(ji,jj) = hbdy(ji,jj)
[9019]550         ENDDO
[11574]551      ENDDO
[6140]552      !
[9019]553      ! --- East --- !
[11574]554      istart = jpiglo - nbghostcells
555      iend   = jpiglo - 1
556      DO ji = mi0(istart), mi1(iend)
[9019]557         DO jj = 1, jpj
[11574]558            ssha(ji,jj) = hbdy(ji,jj)
[9019]559         ENDDO
[11574]560      ENDDO
[6140]561      !
[9019]562      ! --- South --- !
[11574]563      jstart = 2
564      jend   = 1 + nbghostcells
565      DO jj = mj0(jstart), mj1(jend)
566         DO ji = 1, jpi
567            ssha(ji,jj) = hbdy(ji,jj)
[9019]568         ENDDO
[11574]569      ENDDO
[6140]570      !
[9019]571      ! --- North --- !
[11574]572      jstart = jpjglo - nbghostcells
573      jend   = jpjglo - 1
574      DO jj = mj0(jstart), mj1(jend)
575         DO ji = 1, jpi
576            ssha(ji,jj) = hbdy(ji,jj)
[9019]577         ENDDO
[11574]578      ENDDO
[6140]579      !
[2486]580   END SUBROUTINE Agrif_ssh
581
[6140]582
[4486]583   SUBROUTINE Agrif_ssh_ts( jn )
[4292]584      !!----------------------------------------------------------------------
585      !!                  ***  ROUTINE Agrif_ssh_ts  ***
586      !!---------------------------------------------------------------------- 
[4486]587      INTEGER, INTENT(in) ::   jn
[4292]588      !!
[11574]589      INTEGER :: ji, jj
590      INTEGER  :: istart, iend, jstart, jend
[4292]591      !!---------------------------------------------------------------------- 
[9031]592      !
593      IF( Agrif_Root() )   RETURN
594      !
[9116]595      ! --- West --- !
[11574]596      istart = 2
597      iend   = 1+nbghostcells
598      DO ji = mi0(istart), mi1(iend)
[6140]599         DO jj = 1, jpj
[11574]600            ssha_e(ji,jj) = hbdy(ji,jj)
[9116]601         ENDDO
[11574]602      ENDDO
[6140]603      !
[9116]604      ! --- East --- !
[11574]605      istart = jpiglo - nbghostcells
606      iend   = jpiglo - 1
607      DO ji = mi0(istart), mi1(iend)
[6140]608         DO jj = 1, jpj
[11574]609            ssha_e(ji,jj) = hbdy(ji,jj)
[9116]610         ENDDO
[11574]611      ENDDO
[6140]612      !
[9116]613      ! --- South --- !
[11574]614      jstart = 2
615      jend   = 1+nbghostcells
616      DO jj = mj0(jstart), mj1(jend)
617         DO ji = 1, jpi
618            ssha_e(ji,jj) = hbdy(ji,jj)
[9116]619         ENDDO
[11574]620      ENDDO
[6140]621      !
[9116]622      ! --- North --- !
[11574]623      jstart = jpjglo - nbghostcells
624      jend   = jpjglo - 1
625      DO jj = mj0(jstart), mj1(jend)
626         DO ji = 1, jpi
627            ssha_e(ji,jj) = hbdy(ji,jj)
[9116]628         ENDDO
[11574]629      ENDDO
[6140]630      !
[4292]631   END SUBROUTINE Agrif_ssh_ts
632
[9019]633   SUBROUTINE Agrif_avm
[4292]634      !!----------------------------------------------------------------------
[9019]635      !!                  ***  ROUTINE Agrif_avm  ***
[5656]636      !!---------------------------------------------------------------------- 
637      REAL(wp) ::   zalpha
[6140]638      !!---------------------------------------------------------------------- 
[5656]639      !
[9031]640      IF( Agrif_Root() )   RETURN
[6140]641      !
[9031]642      zalpha = 1._wp ! JC: proper time interpolation impossible 
643                     ! => use last available value from parent
644      !
645      Agrif_SpecialValue    = 0.e0
[5656]646      Agrif_UseSpecialValue = .TRUE.
[6140]647      !
[9019]648      CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )       
[6140]649      !
[5656]650      Agrif_UseSpecialValue = .FALSE.
651      !
[9019]652   END SUBROUTINE Agrif_avm
[6140]653   
[5656]654
[11574]655   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
[6140]656      !!----------------------------------------------------------------------
[9019]657      !!                  *** ROUTINE interptsn ***
[6140]658      !!----------------------------------------------------------------------
659      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
660      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
661      LOGICAL                                     , INTENT(in   ) ::   before
[5656]662      !
[11574]663      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
664      INTEGER  ::   N_in, N_out
[9031]665      ! vertical interpolation:
[11741]666      REAL(wp) :: zhtot
[11625]667      REAL(wp), DIMENSION(k1:k2,1:jpts) :: tabin
[9031]668      REAL(wp), DIMENSION(k1:k2) :: h_in
[9116]669      REAL(wp), DIMENSION(1:jpk) :: h_out
[11590]670      !!----------------------------------------------------------------------
[9031]671
[9019]672      IF( before ) THEN         
[9031]673         DO jn = 1,jpts
674            DO jk=k1,k2
675               DO jj=j1,j2
676                 DO ji=i1,i2
677                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
678                 END DO
679              END DO
680           END DO
681        END DO
682
683# if defined key_vertical
[11741]684        ! Interpolate thicknesses
685        ! Warning: these are masked, hence extrapolated prior interpolation.
[11769]686        DO jk=k1,k2
[9031]687           DO jj=j1,j2
688              DO ji=i1,i2
[11741]689                  ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
[9031]690              END DO
691           END DO
692        END DO
[11741]693
694        ! Extrapolate thicknesses in partial bottom cells:
695        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
696        IF (ln_zps) THEN
697           DO jj=j1,j2
698              DO ji=i1,i2
699                  jk = mbkt(ji,jj)
700                  ptab(ji,jj,jk,jpts+1) = 0._wp
701              END DO
702           END DO           
703        END IF
704     
705        ! Save ssh at last level:
706        IF (.NOT.ln_linssh) THEN
707           ptab(i1:i2,j1:j2,k2,jpts+1) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 
708        ELSE
709           ptab(i1:i2,j1:j2,k2,jpts+1) = 0._wp
710        END IF     
[9031]711# endif
712      ELSE 
713
[11741]714# if defined key_vertical 
715         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
716           
[9031]717         DO jj=j1,j2
718            DO ji=i1,i2
[11741]719               tsa(ji,jj,:,:) = 0._wp
720               N_in = mbkt_parent(ji,jj)
721               zhtot = 0._wp
722               DO jk=1,N_in !k2 = jpk of parent grid
723                  IF (jk==N_in) THEN
724                     h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot
725                  ELSE
726                     h_in(jk) = ptab(ji,jj,jk,n2)
727                  ENDIF
728                  zhtot = zhtot + h_in(jk)
[9031]729                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
730               END DO
731               N_out = 0
732               DO jk=1,jpk ! jpk of child grid
[11769]733                  IF (tmask(ji,jj,jk) == 0._wp) EXIT
[9031]734                  N_out = N_out + 1
[11590]735                  h_out(jk) = e3t_a(ji,jj,jk)
[9031]736               ENDDO
[11625]737               IF (N_in*N_out > 0) THEN
[11741]738                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tsa(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
[9031]739               ENDIF
740            ENDDO
741         ENDDO
742# else
[5656]743         !
[9759]744         DO jn=1, jpts
[11741]745            tsa(i1:i2,j1:j2,1:jpk,jn)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
[9759]746         END DO
[11741]747# endif
[9116]748
[5656]749      ENDIF
750      !
751   END SUBROUTINE interptsn
752
[11574]753   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before )
[5656]754      !!----------------------------------------------------------------------
[4292]755      !!                  ***  ROUTINE interpsshn  ***
756      !!---------------------------------------------------------------------- 
[6140]757      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
758      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
759      LOGICAL                         , INTENT(in   ) ::   before
760      !
[5656]761      !!---------------------------------------------------------------------- 
762      !
763      IF( before) THEN
764         ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
765      ELSE
[11574]766         hbdy(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
[5656]767      ENDIF
768      !
769   END SUBROUTINE interpsshn
770
[11590]771   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
[6140]772      !!----------------------------------------------------------------------
[9019]773      !!                  *** ROUTINE interpun ***
[9031]774      !!---------------------------------------------   
775      !!
776      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
777      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
778      LOGICAL, INTENT(in) :: before
779      !!
780      INTEGER :: ji,jj,jk
[11741]781      REAL(wp) :: zrhoy, zhtot
[9031]782      ! vertical interpolation:
783      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
784      REAL(wp), DIMENSION(1:jpk) :: h_out
[11590]785      INTEGER  :: N_in, N_out
[9031]786      REAL(wp) :: h_diff
787      !!---------------------------------------------   
[5656]788      !
[9031]789      IF (before) THEN
790         DO jk=1,jpk
791            DO jj=j1,j2
792               DO ji=i1,i2
793                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 
794# if defined key_vertical
[11741]795                  ! Interpolate thicknesses (masked for subsequent extrapolation)
796                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)
[9031]797# endif
798               END DO
799            END DO
[5656]800         END DO
[11741]801# if defined key_vertical
802         ! Extrapolate thicknesses in partial bottom cells:
803         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
804         IF (ln_zps) THEN
805            DO jj=j1,j2
806               DO ji=i1,i2
807                  jk = mbku(ji,jj)
808                  ptab(ji,jj,jk,2) = 0._wp
809               END DO
810            END DO           
811         END IF
812        ! Save ssh at last level:
813        ptab(i1:i2,j1:j2,k2,2) = 0._wp
814        IF (.NOT.ln_linssh) THEN
815           ! This vertical sum below should be replaced by the sea-level at U-points (optimization):
816           DO jk=1,jpk
817              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3u_n(i1:i2,j1:j2,jk) * umask(i1:i2,j1:j2,jk)
818           END DO
819           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hu_0(i1:i2,j1:j2)
820        END IF 
821# endif
822         !
[5656]823      ELSE
[9031]824         zrhoy = Agrif_rhoy()
825# if defined key_vertical
826! VERTICAL REFINEMENT BEGIN
827
[11741]828         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
829
[9031]830         DO ji=i1,i2
831            DO jj=j1,j2
[11741]832               ua(ji,jj,:) = 0._wp
833               N_in = mbku_parent(ji,jj)
834               zhtot = 0._wp
835               DO jk=1,N_in
836                  IF (jk==N_in) THEN
837                     h_in(jk) = hu0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
838                  ELSE
839                     h_in(jk) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
840                  ENDIF
841                  zhtot = zhtot + h_in(jk)
842                  tabin(jk) = ptab(ji,jj,jk,1)/(e2u(ji,jj)*zrhoy*h_in(jk))
[9031]843              ENDDO
[11741]844                 
[9031]845              N_out = 0
846              DO jk=1,jpk
[11590]847                 if (umask(ji,jj,jk) == 0) EXIT
[9031]848                 N_out = N_out + 1
[11590]849                 h_out(N_out) = e3u_a(ji,jj,jk)
[9031]850              ENDDO
[11741]851              IF (N_in*N_out > 0) THEN
852                 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
[9031]853              ENDIF
854            ENDDO
855         ENDDO
856
857# else
[6140]858         DO jk = 1, jpkm1
[5656]859            DO jj=j1,j2
[9031]860               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) )
[5656]861            END DO
862         END DO
[9031]863# endif
864
[5656]865      ENDIF
866      !
867   END SUBROUTINE interpun
868
[11590]869   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
[6140]870      !!----------------------------------------------------------------------
[9019]871      !!                  *** ROUTINE interpvn ***
[6140]872      !!----------------------------------------------------------------------
[5656]873      !
[9031]874      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
875      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
876      LOGICAL, INTENT(in) :: before
877      !
878      INTEGER :: ji,jj,jk
879      REAL(wp) :: zrhox
880      ! vertical interpolation:
881      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
882      REAL(wp), DIMENSION(1:jpk) :: h_out
[11590]883      INTEGER  :: N_in, N_out
[11741]884      REAL(wp) :: h_diff, zhtot
[9031]885      !!---------------------------------------------   
[5656]886      !     
[9031]887      IF (before) THEN         
888         DO jk=k1,k2
889            DO jj=j1,j2
890               DO ji=i1,i2
891                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
892# if defined key_vertical
[11741]893                  ! Interpolate thicknesses (masked for subsequent extrapolation)
[9031]894                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
895# endif
896               END DO
897            END DO
[5656]898         END DO
[11741]899# if defined key_vertical
900         ! Extrapolate thicknesses in partial bottom cells:
901         ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
902         IF (ln_zps) THEN
903            DO jj=j1,j2
904               DO ji=i1,i2
905                  jk = mbkv(ji,jj)
906                  ptab(ji,jj,jk,2) = 0._wp
907               END DO
908            END DO           
909         END IF
910        ! Save ssh at last level:
911        ptab(i1:i2,j1:j2,k2,2) = 0._wp
912        IF (.NOT.ln_linssh) THEN
913           ! This vertical sum below should be replaced by the sea-level at V-points (optimization):
914           DO jk=1,jpk
915              ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) + e3v_n(i1:i2,j1:j2,jk) * vmask(i1:i2,j1:j2,jk)
916           END DO
917           ptab(i1:i2,j1:j2,k2,2) = ptab(i1:i2,j1:j2,k2,2) - hv_0(i1:i2,j1:j2)
918        END IF 
919# endif
[9031]920      ELSE       
921         zrhox = Agrif_rhox()
922# if defined key_vertical
923
[11741]924         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
925
[9031]926         DO jj=j1,j2
927            DO ji=i1,i2
[11741]928               va(ji,jj,:) = 0._wp
929               N_in = mbkv_parent(ji,jj)
930               zhtot = 0._wp
931               DO jk=1,N_in
932                  IF (jk==N_in) THEN
933                     h_in(jk) = hv0_parent(ji,jj) + ptab(ji,jj,k2,2) - zhtot
934                  ELSE
935                     h_in(jk) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
936                  ENDIF
937                  zhtot = zhtot + h_in(jk)
938                  tabin(jk) = ptab(ji,jj,jk,1)/(e1v(ji,jj)*zrhox*h_in(jk))
939              ENDDO
[9031]940         
941               N_out = 0
942               DO jk=1,jpk
[11590]943                  if (vmask(ji,jj,jk) == 0) EXIT
[9031]944                  N_out = N_out + 1
[11590]945                  h_out(N_out) = e3v_a(ji,jj,jk)
[9031]946               END DO
[11741]947               IF (N_in*N_out > 0) THEN
948                  call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
[9031]949               ENDIF
950            END DO
951         END DO
952# else
[6140]953         DO jk = 1, jpkm1
[9031]954            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) )
[5656]955         END DO
[9031]956# endif
[5656]957      ENDIF
958      !       
959   END SUBROUTINE interpvn
[636]960
[11574]961   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before)
[1605]962      !!----------------------------------------------------------------------
[5656]963      !!                  ***  ROUTINE interpunb  ***
[1605]964      !!---------------------------------------------------------------------- 
[6140]965      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
966      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
967      LOGICAL                         , INTENT(in   ) ::   before
968      !
969      INTEGER  ::   ji, jj
970      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
[1605]971      !!---------------------------------------------------------------------- 
[5656]972      !
[6140]973      IF( before ) THEN
974         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
[5656]975      ELSE
976         zrhoy = Agrif_Rhoy()
977         zrhot = Agrif_rhot()
978         ! Time indexes bounds for integration
979         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
980         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot     
[11574]981         !
982         DO ji = i1, i2
983            DO jj = j1, j2
[11603]984               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
985                  IF    ( utint_stage(ji,jj) == 1  ) THEN
986                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
987                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
988                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN
989                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
990                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
991                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN               
992                     ztcoeff = 1._wp
993                  ELSE
994                     ztcoeff = 0._wp
995                  ENDIF
996                  !   
997                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
998                  !           
999                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
1000                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
1001                  ENDIF
1002                  !
[11574]1003                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
1004               ENDIF
1005            END DO
1006         END DO
1007      END IF
[5656]1008      !
1009   END SUBROUTINE interpunb
[636]1010
[6140]1011
[11574]1012   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before )
[1605]1013      !!----------------------------------------------------------------------
[5656]1014      !!                  ***  ROUTINE interpvnb  ***
[1605]1015      !!---------------------------------------------------------------------- 
[6140]1016      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1017      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1018      LOGICAL                         , INTENT(in   ) ::   before
1019      !
[11574]1020      INTEGER  ::   ji, jj
[6140]1021      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
[1605]1022      !!---------------------------------------------------------------------- 
[5656]1023      !
[6140]1024      IF( before ) THEN
1025         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
[5656]1026      ELSE
1027         zrhox = Agrif_Rhox()
1028         zrhot = Agrif_rhot()
1029         ! Time indexes bounds for integration
1030         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
[11574]1031         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
1032         !     
1033         DO ji = i1, i2
1034            DO jj = j1, j2
[11603]1035               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
1036                  IF    ( vtint_stage(ji,jj) == 1  ) THEN
1037                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
1038                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
1039                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
1040                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
1041                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
1042                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN               
1043                     ztcoeff = 1._wp
1044                  ELSE
1045                     ztcoeff = 0._wp
1046                  ENDIF
1047                  !   
1048                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
1049                  !           
1050                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
1051                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
1052                  ENDIF
1053                  !
[11574]1054                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
1055               ENDIF
1056            END DO
1057         END DO         
[5656]1058      ENDIF
1059      !
1060   END SUBROUTINE interpvnb
[390]1061
[6140]1062
[11574]1063   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before )
[1605]1064      !!----------------------------------------------------------------------
[5656]1065      !!                  ***  ROUTINE interpub2b  ***
[1605]1066      !!---------------------------------------------------------------------- 
[6140]1067      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1068      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1069      LOGICAL                         , INTENT(in   ) ::   before
1070      !
1071      INTEGER  ::   ji,jj
[11574]1072      REAL(wp) ::   zrhot, zt0, zt1, zat
[1605]1073      !!---------------------------------------------------------------------- 
[5656]1074      IF( before ) THEN
[9031]1075         IF ( ln_bt_fw ) THEN
1076            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1077         ELSE
1078            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1079         ENDIF
[5656]1080      ELSE
1081         zrhot = Agrif_rhot()
1082         ! Time indexes bounds for integration
1083         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1084         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1085         ! Polynomial interpolation coefficients:
[6140]1086         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1087            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
[11574]1088         !
1089         ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
1090         !
1091         ! Update interpolation stage:
1092         utint_stage(i1:i2,j1:j2) = 1
[5656]1093      ENDIF
1094      !
1095   END SUBROUTINE interpub2b
[6140]1096   
[636]1097
[11574]1098   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before )
[4292]1099      !!----------------------------------------------------------------------
[5656]1100      !!                  ***  ROUTINE interpvb2b  ***
[4292]1101      !!---------------------------------------------------------------------- 
[6140]1102      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1103      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1104      LOGICAL                         , INTENT(in   ) ::   before
1105      !
1106      INTEGER ::   ji,jj
[11574]1107      REAL(wp) ::   zrhot, zt0, zt1, zat
[4292]1108      !!---------------------------------------------------------------------- 
[5656]1109      !
1110      IF( before ) THEN
[9031]1111         IF ( ln_bt_fw ) THEN
1112            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1113         ELSE
1114            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1115         ENDIF
[5656]1116      ELSE     
1117         zrhot = Agrif_rhot()
1118         ! Time indexes bounds for integration
1119         zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
1120         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1121         ! Polynomial interpolation coefficients:
[6140]1122         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
1123            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
[5656]1124         !
[11574]1125         vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
1126         !
1127         ! update interpolation stage:
1128         vtint_stage(i1:i2,j1:j2) = 1
[5656]1129      ENDIF
1130      !     
1131   END SUBROUTINE interpvb2b
[4292]1132
[6140]1133
[12119]1134   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before )
[5656]1135      !!----------------------------------------------------------------------
1136      !!                  ***  ROUTINE interpe3t  ***
1137      !!---------------------------------------------------------------------- 
[6140]1138      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
[5656]1139      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
[6140]1140      LOGICAL                              , INTENT(in   ) :: before
[5656]1141      !
1142      INTEGER :: ji, jj, jk
1143      !!---------------------------------------------------------------------- 
1144      !   
[6140]1145      IF( before ) THEN
1146         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
[5656]1147      ELSE
[9019]1148         !
[6140]1149         DO jk = k1, k2
1150            DO jj = j1, j2
1151               DO ji = i1, i2
[9019]1152                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
[12119]1153                     WRITE(numout,*) ' Agrif error for e3t_0: parent , child, i, j, k ',  & 
1154                     &                 ptab(ji,jj,jk), tmask(ji,jj,jk) * e3t_0(ji,jj,jk), &
1155                     &                 ji+nimpp-1, jj+njmpp-1, jk
1156                     kindic_agr = kindic_agr + 1
[5656]1157                  ENDIF
1158               END DO
1159            END DO
1160         END DO
[6140]1161         !
[5656]1162      ENDIF
1163      !
1164   END SUBROUTINE interpe3t
1165
[6140]1166
[9031]1167   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
[4486]1168      !!----------------------------------------------------------------------
[5656]1169      !!                  ***  ROUTINE interavm  ***
[4486]1170      !!---------------------------------------------------------------------- 
[9116]1171      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2
[9031]1172      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab
[9116]1173      LOGICAL                                    , INTENT(in   ) ::   before
[11802]1174      !
1175      INTEGER  :: ji, jj, jk
1176      INTEGER  :: N_in, N_out
1177      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in
1178      REAL(wp), DIMENSION(1:jpk) :: z_out
[4486]1179      !!---------------------------------------------------------------------- 
[5656]1180      !     
[9031]1181      IF (before) THEN         
1182         DO jk=k1,k2
1183            DO jj=j1,j2
1184              DO ji=i1,i2
1185                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1186              END DO
1187           END DO
1188        END DO
[11802]1189
1190# if defined key_vertical
1191        ! Interpolate thicknesses
1192        ! Warning: these are masked, hence extrapolated prior interpolation.
[9031]1193        DO jk=k1,k2
1194           DO jj=j1,j2
1195              DO ji=i1,i2
[11802]1196                  ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
[9031]1197              END DO
1198           END DO
1199        END DO
[11802]1200
1201        ! Extrapolate thicknesses in partial bottom cells:
1202        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
1203        IF (ln_zps) THEN
1204           DO jj=j1,j2
1205              DO ji=i1,i2
1206                  jk = mbkt(ji,jj)
1207                  ptab(ji,jj,jk,2) = 0._wp
1208              END DO
1209           END DO           
1210        END IF
1211     
1212        ! Save ssh at last level:
1213        IF (.NOT.ln_linssh) THEN
1214           ptab(i1:i2,j1:j2,k2,2) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 
1215        ELSE
1216           ptab(i1:i2,j1:j2,k2,2) = 0._wp
1217        END IF     
1218# endif
[9031]1219      ELSE 
1220#ifdef key_vertical         
[11802]1221         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp 
1222         avm_k(i1:i2,j1:j2,k1:k2) = 0._wp
1223           
1224         DO jj = j1, j2
1225            DO ji =i1, i2
1226               N_in = mbkt_parent(ji,jj)
1227               IF ( tmask(ji,jj,1) == 0._wp) N_in = 0
1228               z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2)
1229               DO jk = N_in, 1, -1  ! Parent vertical grid               
1230                     z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2)
1231                    tabin(jk) = ptab(ji,jj,jk,1)
[9031]1232               END DO
[12152]1233               N_out = mbkt(ji,jj) 
1234               DO jk = 1, N_out        ! Child vertical grid
1235                  z_out(jk) = gdepw_n(ji,jj,jk)
[9031]1236               ENDDO
[11802]1237               IF (N_in*N_out > 0) THEN
1238                  CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1)
[9031]1239               ENDIF
1240            ENDDO
1241         ENDDO
1242#else
1243         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1244#endif
[5656]1245      ENDIF
1246      !
1247   END SUBROUTINE interpavm
[4486]1248
[11741]1249# if defined key_vertical
1250   SUBROUTINE interpmbkt( ptab, i1, i2, j1, j2, before )
1251      !!----------------------------------------------------------------------
1252      !!                  ***  ROUTINE interpsshn  ***
1253      !!---------------------------------------------------------------------- 
1254      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1255      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1256      LOGICAL                         , INTENT(in   ) ::   before
1257      !
1258      !!---------------------------------------------------------------------- 
1259      !
1260      IF( before) THEN
1261         ptab(i1:i2,j1:j2) = REAL(mbkt(i1:i2,j1:j2),wp)
1262      ELSE
1263         mbkt_parent(i1:i2,j1:j2) = NINT(ptab(i1:i2,j1:j2))
1264      ENDIF
1265      !
1266   END SUBROUTINE interpmbkt
1267
1268   SUBROUTINE interpht0( ptab, i1, i2, j1, j2, before )
1269      !!----------------------------------------------------------------------
1270      !!                  ***  ROUTINE interpsshn  ***
1271      !!---------------------------------------------------------------------- 
1272      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1273      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
1274      LOGICAL                         , INTENT(in   ) ::   before
1275      !
1276      !!---------------------------------------------------------------------- 
1277      !
1278      IF( before) THEN
1279         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2)
1280      ELSE
1281         ht0_parent(i1:i2,j1:j2) = ptab(i1:i2,j1:j2)
1282      ENDIF
1283      !
1284   END SUBROUTINE interpht0
1285#endif
1286
[390]1287#else
[1605]1288   !!----------------------------------------------------------------------
1289   !!   Empty module                                          no AGRIF zoom
1290   !!----------------------------------------------------------------------
[636]1291CONTAINS
[9570]1292   SUBROUTINE Agrif_OCE_Interp_empty
1293      WRITE(*,*)  'agrif_oce_interp : You should not have seen this print! error?'
1294   END SUBROUTINE Agrif_OCE_Interp_empty
[390]1295#endif
[1605]1296
1297   !!======================================================================
[9570]1298END MODULE agrif_oce_interp
Note: See TracBrowser for help on using the repository browser.