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_update.F90 in NEMO/releases/r4.0/r4.0-HEAD/src/NST – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/NST/agrif_oce_update.F90 @ 15564

Last change on this file since 15564 was 15564, checked in by jchanut, 3 years ago

#2751, fixes small constancy preservation issue with AGRIF

  • Property svn:keywords set to Id
File size: 55.3 KB
Line 
1#define TWO_WAY        /* TWO WAY NESTING */
2#undef DECAL_FEEDBACK  /* SEPARATION of INTERFACES*/
3#undef VOL_REFLUX      /* VOLUME REFLUXING*/
4 
5MODULE agrif_oce_update
6   !!======================================================================
7   !!                   ***  MODULE  agrif_oce_interp  ***
8   !! AGRIF: update package for the ocean dynamics (OPA)
9   !!======================================================================
10   !! History :  2.0  !  2002-06  (L. Debreu)  Original code
11   !!            3.2  !  2009-04  (R. Benshila)
12   !!            3.6  !  2014-09  (R. Benshila)
13   !!----------------------------------------------------------------------
14#if defined key_agrif 
15   !!----------------------------------------------------------------------
16   !!   'key_agrif'                                              AGRIF zoom
17   !!----------------------------------------------------------------------
18   USE par_oce
19   USE oce
20   USE dom_oce
21   USE zdf_oce        ! vertical physics: ocean variables
22   USE agrif_oce
23   !
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE domvvl         ! Need interpolation routines
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh
32   PUBLIC   Update_Scales
33
34   !!----------------------------------------------------------------------
35   !! NEMO/NST 4.0 , NEMO Consortium (2018)
36   !! $Id$
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE Agrif_Update_Tra( )
42      !!----------------------------------------------------------------------
43      !!                   *** ROUTINE Agrif_Update_Tra ***
44      !!----------------------------------------------------------------------
45      !
46      IF (Agrif_Root()) RETURN
47      !
48#if defined TWO_WAY 
49      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed()
50
51      Agrif_UseSpecialValueInUpdate = .TRUE.
52      Agrif_SpecialValueFineGrid    = 0._wp
53      !
54# if ! defined DECAL_FEEDBACK
55      CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
56! near boundary update:
57!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
58# else
59      CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
60! near boundary update:
61!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
62# endif
63      !
64      Agrif_UseSpecialValueInUpdate = .FALSE.
65      !
66#endif
67      !
68   END SUBROUTINE Agrif_Update_Tra
69
70   SUBROUTINE Agrif_Update_Dyn( )
71      !!----------------------------------------------------------------------
72      !!                   *** ROUTINE Agrif_Update_Dyn ***
73      !!----------------------------------------------------------------------
74      !
75      IF (Agrif_Root()) RETURN
76      !
77#if defined TWO_WAY
78      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed()
79
80      Agrif_UseSpecialValueInUpdate = .FALSE.
81      Agrif_SpecialValueFineGrid = 0.
82      !     
83# if ! defined DECAL_FEEDBACK
84      CALL Agrif_Update_Variable(un_update_id,procname = updateU)
85      CALL Agrif_Update_Variable(vn_update_id,procname = updateV)
86! near boundary update:
87!      CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
88!      CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)
89# else
90      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
91      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
92! near boundary update:
93!      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
94!      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
95# endif
96
97# if ! defined DECAL_FEEDBACK
98      CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
99      CALL Agrif_Update_Variable(e2v_id,procname = updateV2d) 
100# else
101      CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
102      CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d) 
103# endif
104      !
105# if ! defined DECAL_FEEDBACK
106      ! Account for updated thicknesses at boundary edges
107      IF (.NOT.ln_linssh) THEN
108!         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)
109!         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)
110      ENDIF
111# endif
112      !
113      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN
114         ! Update time integrated transports
115#  if ! defined DECAL_FEEDBACK
116         CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
117         CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
118#  else
119         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
120         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
121#  endif
122      END IF
123#endif
124      !
125   END SUBROUTINE Agrif_Update_Dyn
126
127   SUBROUTINE Agrif_Update_ssh( )
128      !!---------------------------------------------
129      !!   *** ROUTINE Agrif_Update_ssh ***
130      !!---------------------------------------------
131      !
132      IF (Agrif_Root()) RETURN
133      !
134#if defined TWO_WAY
135      !
136      Agrif_UseSpecialValueInUpdate = .TRUE.
137      Agrif_SpecialValueFineGrid = 0.
138# if ! defined DECAL_FEEDBACK
139      CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
140# else
141      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
142# endif
143      !
144      Agrif_UseSpecialValueInUpdate = .FALSE.
145      !
146#  if defined VOL_REFLUX
147      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
148         ! Refluxing on ssh:
149#  if defined DECAL_FEEDBACK
150         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)
151         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv)
152#  else
153         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu)
154         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv)
155#  endif
156      END IF
157#  endif
158      !
159#endif
160      !
161   END SUBROUTINE Agrif_Update_ssh
162
163
164   SUBROUTINE Agrif_Update_Tke( )
165      !!---------------------------------------------
166      !!   *** ROUTINE Agrif_Update_Tke ***
167      !!---------------------------------------------
168      !!
169      !
170      IF (Agrif_Root()) RETURN
171      !       
172#  if defined TWO_WAY
173
174      Agrif_UseSpecialValueInUpdate = .TRUE.
175      Agrif_SpecialValueFineGrid = 0.
176
177      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  )
178      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
179      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
180
181      Agrif_UseSpecialValueInUpdate = .FALSE.
182
183#  endif
184     
185   END SUBROUTINE Agrif_Update_Tke
186
187
188   SUBROUTINE Agrif_Update_vvl( )
189      !!---------------------------------------------
190      !!   *** ROUTINE Agrif_Update_vvl ***
191      !!---------------------------------------------
192      !
193      IF (Agrif_Root()) RETURN
194      !
195#if defined TWO_WAY 
196      !
197      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step()
198      !
199      Agrif_UseSpecialValueInUpdate = .TRUE.
200      Agrif_SpecialValueFineGrid = 0.
201      !
202      ! No interface separation here, update vertical grid at T points
203      ! everywhere over the overlapping regions (one account for refluxing in that case):
204      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t) 
205      !
206      Agrif_UseSpecialValueInUpdate = .FALSE.
207      !
208      CALL Agrif_ChildGrid_To_ParentGrid()
209      CALL dom_vvl_update_UVF
210      CALL Agrif_ParentGrid_To_ChildGrid()
211      !
212#endif
213      !
214   END SUBROUTINE Agrif_Update_vvl
215
216   SUBROUTINE dom_vvl_update_UVF
217      !!---------------------------------------------
218      !!       *** ROUTINE dom_vvl_update_UVF ***
219      !!---------------------------------------------
220      !!
221      INTEGER :: jk
222      REAL(wp):: zcoef
223      !!---------------------------------------------
224
225      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', &
226                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step()
227
228      ! Save "old" scale factor (prior update) for subsequent asselin correction
229      ! of prognostic variables
230      ! -----------------------
231      !
232      e3u_a(:,:,:) = e3u_n(:,:,:)
233      e3v_a(:,:,:) = e3v_n(:,:,:)
234!      ua(:,:,:) = e3u_b(:,:,:)
235!      va(:,:,:) = e3v_b(:,:,:)
236      hu_a(:,:) = hu_n(:,:)
237      hv_a(:,:) = hv_n(:,:)
238
239      ! 1) NOW fields
240      !--------------
241     
242         ! Vertical scale factor interpolations
243         ! ------------------------------------
244      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' )
245      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' )
246      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' )
247
248      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
249      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
250
251         ! Update total depths:
252         ! --------------------
253      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points
254      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points
255      DO jk = 1, jpkm1
256         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
257         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
258      END DO
259      !                                        ! Inverse of the local depth
260      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
261      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
262
263
264      ! 2) BEFORE fields:
265      !------------------
266      IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
267         !
268         ! Vertical scale factor interpolations
269         ! ------------------------------------
270         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  )
271         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  )
272
273         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
274         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
275
276         ! Update total depths:
277         ! --------------------
278         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points
279         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points
280         DO jk = 1, jpkm1
281            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
282            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
283         END DO
284         !                                     ! Inverse of the local depth
285         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
286         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
287      ENDIF
288      !
289   END SUBROUTINE dom_vvl_update_UVF
290
291#if defined key_vertical
292
293   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
294      !!----------------------------------------------------------------------
295      !!           *** ROUTINE updateT ***
296      !!---------------------------------------------
297      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
298      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
299      LOGICAL, INTENT(in) :: before
300      !!
301      INTEGER :: ji,jj,jk,jn
302      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
303      REAL(wp) :: h_in(k1:k2)
304      REAL(wp) :: h_out(1:jpk)
305      INTEGER  :: N_in, N_out
306      REAL(wp) :: zrho_xy, h_diff
307      REAL(wp) :: tabin(k1:k2,n1:n2)
308      !!---------------------------------------------
309      !
310      IF (before) THEN
311         AGRIF_SpecialValue = -999._wp
312         zrho_xy = Agrif_rhox() * Agrif_rhoy() 
313         DO jn = n1,n2-1
314            DO jk=k1,k2
315               DO jj=j1,j2
316                  DO ji=i1,i2
317                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) &
318                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
319                  END DO
320               END DO
321            END DO
322         END DO
323         DO jk=k1,k2
324            DO jj=j1,j2
325               DO ji=i1,i2
326                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
327                                           + (tmask(ji,jj,jk)-1)*999._wp
328               END DO
329            END DO
330         END DO
331      ELSE
332         tabres_child(:,:,:,:) = 0.
333         AGRIF_SpecialValue = 0._wp
334         DO jj=j1,j2
335            DO ji=i1,i2
336               N_in = 0
337               DO jk=k1,k2 !k2 = jpk of child grid
338                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT
339                  N_in = N_in + 1
340                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
341                  h_in(N_in) = tabres(ji,jj,jk,n2)
342               ENDDO
343               N_out = 0
344               DO jk=1,jpk ! jpk of parent grid
345                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
346                  N_out = N_out + 1
347                  h_out(N_out) = e3t_n(ji,jj,jk) 
348               ENDDO
349               IF (N_in > 0) THEN !Remove this?
350                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
351                  IF (h_diff < -1.e-4) THEN
352                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
353                     print *,h_in(1:N_in)
354                     print *,h_out(1:N_out)
355                     STOP
356                  ENDIF
357                  DO jn=n1,n2-1
358                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
359                  ENDDO
360               ENDIF
361            ENDDO
362         ENDDO
363
364         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
365            ! Add asselin part
366            DO jn = n1,n2-1
367               DO jk=1,jpk
368                  DO jj=j1,j2
369                     DO ji=i1,i2
370                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
371                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
372                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
373                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
374                        ENDIF
375                     ENDDO
376                  ENDDO
377               ENDDO
378            ENDDO
379         ENDIF
380         DO jn = n1,n2-1
381            DO jk=1,jpk
382               DO jj=j1,j2
383                  DO ji=i1,i2
384                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
385                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
386                     END IF
387                  END DO
388               END DO
389            END DO
390         END DO
391      ENDIF
392      !
393   END SUBROUTINE updateTS
394
395# else
396
397   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
398      !!---------------------------------------------
399      !!           *** ROUTINE updateT ***
400      !!---------------------------------------------
401      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
402      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
403      LOGICAL, INTENT(in) :: before
404      !!
405      INTEGER :: ji,jj,jk,jn
406      REAL(wp) :: ze3b, ztb, ztnu, ztno
407      !!---------------------------------------------
408      !
409      IF (before) THEN
410         DO jn = 1,jpts
411            DO jk=k1,k2
412               DO jj=j1,j2
413                  DO ji=i1,i2
414!> jc tmp
415                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
416!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
417!< jc tmp
418                  END DO
419               END DO
420            END DO
421         END DO
422      ELSE
423!> jc tmp
424         DO jn = 1,jpts
425            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
426                                         & * tmask(i1:i2,j1:j2,k1:k2)
427         ENDDO
428!< jc tmp
429         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
430            ! Add asselin part
431            DO jn = 1,jpts
432               DO jk = k1, k2
433                  DO jj = j1, j2
434                     DO ji = i1, i2
435                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
436                           ze3b = e3t_b(ji,jj,jk) & ! Recover e3t_b before update
437                                & - rn_atfp * ( e3t_n(ji,jj,jk) - e3t_a(ji,jj,jk) )
438                           ztb  = tsb(ji,jj,jk,jn) * ze3b 
439                           ztnu = tabres(ji,jj,jk,jn)
440                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
441                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
442                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
443                        ENDIF
444                     END DO
445                  END DO
446               END DO
447            END DO
448         ENDIF
449         DO jn = 1,jpts
450            DO jk=k1,k2
451               DO jj=j1,j2
452                  DO ji=i1,i2
453                     IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
454                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
455                     END IF
456                  END DO
457               END DO
458            END DO
459         END DO
460         !
461         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
462            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts)
463         ENDIF
464         !
465      ENDIF
466      !
467   END SUBROUTINE updateTS
468
469# endif
470
471# if defined key_vertical
472
473   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
474      !!---------------------------------------------
475      !!           *** ROUTINE updateu ***
476      !!---------------------------------------------
477      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
478      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
479      LOGICAL                                     , INTENT(in   ) :: before
480      !
481      INTEGER ::   ji, jj, jk
482      REAL(wp)::   zrhoy
483! VERTICAL REFINEMENT BEGIN
484      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
485      REAL(wp) :: h_in(k1:k2)
486      REAL(wp) :: h_out(1:jpk)
487      INTEGER  :: N_in, N_out
488      REAL(wp) :: h_diff, excess, thick
489      REAL(wp) :: tabin(k1:k2)
490! VERTICAL REFINEMENT END
491      !!---------------------------------------------
492      !
493      IF( before ) THEN
494         zrhoy = Agrif_Rhoy()
495         AGRIF_SpecialValue = -999._wp
496         DO jk=k1,k2
497            DO jj=j1,j2
498               DO ji=i1,i2
499                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
500                                       + (umask(ji,jj,jk)-1)*999._wp
501                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
502                                       + (umask(ji,jj,jk)-1)*999._wp
503               END DO
504            END DO
505         END DO
506      ELSE
507         tabres_child(:,:,:) = 0.
508         AGRIF_SpecialValue = 0._wp
509         DO jj=j1,j2
510            DO ji=i1,i2
511               N_in = 0
512               h_in(:) = 0._wp
513               tabin(:) = 0._wp
514               DO jk=k1,k2 !k2=jpk of child grid
515                  IF( tabres(ji,jj,jk,2) < -900) EXIT
516                  N_in = N_in + 1
517                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
518                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
519               ENDDO
520               N_out = 0
521               DO jk=1,jpk
522                  IF (umask(ji,jj,jk) == 0) EXIT
523                  N_out = N_out + 1
524                  h_out(N_out) = e3u_n(ji,jj,jk)
525               ENDDO
526               IF (N_in * N_out > 0) THEN
527                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
528                  IF (h_diff < -1.e-4) THEN
529!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
530!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
531                     excess = 0._wp
532                     DO jk=N_in,1,-1
533                        thick = MIN(-1*h_diff, h_in(jk))
534                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
535                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
536                        h_diff = h_diff + thick
537                        IF ( h_diff == 0) THEN
538                           N_in = jk
539                           h_in(jk) = h_in(jk) - thick
540                           EXIT
541                        ENDIF
542                     ENDDO
543                  ENDIF
544                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
545                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
546               ENDIF
547            ENDDO
548         ENDDO
549
550         DO jk=1,jpk
551            DO jj=j1,j2
552               DO ji=i1,i2
553                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
554                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
555                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
556                  ENDIF
557                  !
558                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
559               END DO
560            END DO
561         END DO
562      ENDIF
563      !
564   END SUBROUTINE updateu
565
566#else
567
568   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
569      !!---------------------------------------------
570      !!           *** ROUTINE updateu ***
571      !!---------------------------------------------
572      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
573      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
574      LOGICAL                                     , INTENT(in   ) :: before
575      !
576      INTEGER  :: ji, jj, jk
577      REAL(wp) :: zrhoy, ze3b, zub, zunu, zuno
578      !!---------------------------------------------
579      !
580      IF( before ) THEN
581         zrhoy = Agrif_Rhoy()
582         DO jk = k1, k2
583            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
584         END DO
585      ELSE
586         DO jk=k1,k2
587            DO jj=j1,j2
588               DO ji=i1,i2
589                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 
590                  !
591                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
592                     ze3b = e3u_b(ji,jj,jk) & ! Recover e3u_b before update
593                          & - rn_atfp * ( e3u_n(ji,jj,jk) - e3u_a(ji,jj,jk) )
594                     zub  = ub(ji,jj,jk) * ze3b 
595                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
596                     zunu = tabres(ji,jj,jk,1)
597                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &     
598                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
599                  ENDIF
600                  !
601                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
602               END DO
603            END DO
604         END DO
605         !
606         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
607            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
608         ENDIF
609         !
610      ENDIF
611      !
612   END SUBROUTINE updateu
613
614# endif
615
616   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
617      !!---------------------------------------------
618      !!           *** ROUTINE correct_u_bdy ***
619      !!---------------------------------------------
620      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
621      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
622      LOGICAL                                     , INTENT(in   ) :: before
623      INTEGER                                     , INTENT(in)    :: nb, ndir
624      !!
625      LOGICAL :: western_side, eastern_side 
626      !
627      INTEGER  :: jj, jk
628      REAL(wp) :: zcor
629      !!---------------------------------------------
630      !
631      IF( .NOT.before ) THEN
632         !
633         western_side  = (nb == 1).AND.(ndir == 1)
634         eastern_side  = (nb == 1).AND.(ndir == 2)
635         !
636         IF (western_side) THEN
637            DO jj=j1,j2
638               zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
639               un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
640               DO jk=1,jpkm1
641                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)
642               END DO
643            END DO
644         ENDIF
645         !
646         IF (eastern_side) THEN
647            DO jj=j1,j2
648               zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
649               un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
650               DO jk=1,jpkm1
651                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)
652               END DO
653            END DO
654         ENDIF
655         !
656      ENDIF
657      !
658   END SUBROUTINE correct_u_bdy
659
660# if  defined key_vertical
661
662   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
663      !!---------------------------------------------
664      !!           *** ROUTINE updatev ***
665      !!---------------------------------------------
666      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
667      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
668      LOGICAL                                     , INTENT(in   ) :: before
669      !
670      INTEGER  ::   ji, jj, jk
671      REAL(wp) ::   zrhox
672! VERTICAL REFINEMENT BEGIN
673      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
674      REAL(wp) :: h_in(k1:k2)
675      REAL(wp) :: h_out(1:jpk)
676      INTEGER :: N_in, N_out
677      REAL(wp) :: h_diff, excess, thick
678      REAL(wp) :: tabin(k1:k2)
679! VERTICAL REFINEMENT END
680      !!---------------------------------------------     
681      !
682      IF( before ) THEN
683         zrhox = Agrif_Rhox()
684         AGRIF_SpecialValue = -999._wp
685         DO jk=k1,k2
686            DO jj=j1,j2
687               DO ji=i1,i2
688                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
689                                       + (vmask(ji,jj,jk)-1)*999._wp
690                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
691                                       + (vmask(ji,jj,jk)-1)*999._wp
692               END DO
693            END DO
694         END DO
695      ELSE
696         tabres_child(:,:,:) = 0.
697         AGRIF_SpecialValue = 0._wp
698         DO jj=j1,j2
699            DO ji=i1,i2
700               N_in = 0
701               DO jk=k1,k2
702                  IF (tabres(ji,jj,jk,2) < -900) EXIT
703                  N_in = N_in + 1
704                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
705                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
706               ENDDO
707               N_out = 0
708               DO jk=1,jpk
709                  IF (vmask(ji,jj,jk) == 0) EXIT
710                  N_out = N_out + 1
711                  h_out(N_out) = e3v_n(ji,jj,jk)
712               ENDDO
713               IF (N_in * N_out > 0) THEN
714                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
715                  IF (h_diff < -1.e-4) then
716!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
717!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
718                     excess = 0._wp
719                     DO jk=N_in,1,-1
720                        thick = MIN(-1*h_diff, h_in(jk))
721                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
722                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
723                        h_diff = h_diff + thick
724                        IF ( h_diff == 0) THEN
725                           N_in = jk
726                           h_in(jk) = h_in(jk) - thick
727                           EXIT
728                        ENDIF
729                     ENDDO
730                  ENDIF
731                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
732                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
733               ENDIF
734            ENDDO
735         ENDDO
736
737         DO jk=1,jpk
738            DO jj=j1,j2
739               DO ji=i1,i2
740                  !
741                  IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part
742                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
743                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
744                  ENDIF
745                  !
746                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
747               END DO
748            END DO
749         END DO
750      ENDIF
751      !
752   END SUBROUTINE updatev
753
754# else
755
756   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
757      !!---------------------------------------------
758      !!           *** ROUTINE updatev ***
759      !!---------------------------------------------
760      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
761      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
762      LOGICAL                                     , INTENT(in   ) :: before
763      !
764      INTEGER  :: ji, jj, jk
765      REAL(wp) :: zrhox, ze3b, zvb, zvnu, zvno
766      !!---------------------------------------------     
767      !
768      IF (before) THEN
769         zrhox = Agrif_Rhox()
770         DO jk=k1,k2
771            DO jj=j1,j2
772               DO ji=i1,i2
773                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
774               END DO
775            END DO
776         END DO
777      ELSE
778         DO jk=k1,k2
779            DO jj=j1,j2
780               DO ji=i1,i2
781                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
782                  !
783                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
784                     ze3b = e3v_b(ji,jj,jk) & ! Recover e3v_b before update
785                          & - rn_atfp * ( e3v_n(ji,jj,jk) - e3v_a(ji,jj,jk) )
786                     zvb  = vb(ji,jj,jk) * ze3b 
787                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
788                     zvnu = tabres(ji,jj,jk,1)
789                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &     
790                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
791                  ENDIF
792                  !
793                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
794               END DO
795            END DO
796         END DO
797         !
798         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
799            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
800         ENDIF
801         !
802      ENDIF
803      !
804   END SUBROUTINE updatev
805
806# endif
807
808   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
809      !!---------------------------------------------
810      !!           *** ROUTINE correct_u_bdy ***
811      !!---------------------------------------------
812      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
813      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
814      LOGICAL                                     , INTENT(in   ) :: before
815      INTEGER                                     , INTENT(in)    :: nb, ndir
816      !!
817      LOGICAL :: southern_side, northern_side 
818      !
819      INTEGER  :: ji, jk
820      REAL(wp) :: zcor
821      !!---------------------------------------------
822      !
823      IF( .NOT.before ) THEN
824         !
825         southern_side = (nb == 2).AND.(ndir == 1)
826         northern_side = (nb == 2).AND.(ndir == 2)
827         !
828         IF (southern_side) THEN
829            DO ji=i1,i2
830               zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
831               vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
832               DO jk=1,jpkm1
833                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)
834               END DO
835            END DO
836         ENDIF
837         !
838         IF (northern_side) THEN
839            DO ji=i1,i2
840               zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
841               vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
842               DO jk=1,jpkm1
843                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)
844               END DO
845            END DO
846         ENDIF
847         !
848      ENDIF
849      !
850   END SUBROUTINE correct_v_bdy
851
852
853   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
854      !!----------------------------------------------------------------------
855      !!                      *** ROUTINE updateu2d ***
856      !!----------------------------------------------------------------------
857      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
858      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
859      LOGICAL                         , INTENT(in   ) ::   before
860      !!
861      INTEGER  :: ji, jj, jk
862      REAL(wp) :: zrhoy
863      REAL(wp) :: zcorr
864      !!---------------------------------------------
865      !
866      IF( before ) THEN
867         zrhoy = Agrif_Rhoy()
868         DO jj=j1,j2
869            DO ji=i1,i2
870               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
871            END DO
872         END DO
873      ELSE
874         DO jj=j1,j2
875            DO ji=i1,i2
876               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj) 
877               !   
878               ! Update "now" 3d velocities:
879               spgu(ji,jj) = 0._wp
880               DO jk=1,jpkm1
881                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
882               END DO
883               !
884               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
885               DO jk=1,jpkm1             
886                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
887               END DO
888               !
889               ! Update barotropic velocities:
890               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
891                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
892                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
893                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
894                  END IF
895               ENDIF   
896               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
897               !       
898               ! Correct "before" velocities to hold correct bt component:
899               spgu(ji,jj) = 0.e0
900               DO jk=1,jpkm1
901                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
902               END DO
903               !
904               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
905               DO jk=1,jpkm1             
906                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
907               END DO
908               !
909            END DO
910         END DO
911         !
912         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
913            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
914         ENDIF
915      ENDIF
916      !
917   END SUBROUTINE updateu2d
918
919
920   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
921      !!----------------------------------------------------------------------
922      !!                   *** ROUTINE updatev2d ***
923      !!----------------------------------------------------------------------
924      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
925      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
926      LOGICAL                         , INTENT(in   ) ::   before
927      !
928      INTEGER  :: ji, jj, jk
929      REAL(wp) :: zrhox, zcorr
930      !!----------------------------------------------------------------------
931      !
932      IF( before ) THEN
933         zrhox = Agrif_Rhox()
934         DO jj=j1,j2
935            DO ji=i1,i2
936               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
937            END DO
938         END DO
939      ELSE
940         DO jj=j1,j2
941            DO ji=i1,i2
942               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj) 
943               !   
944               ! Update "now" 3d velocities:
945               spgv(ji,jj) = 0.e0
946               DO jk=1,jpkm1
947                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
948               END DO
949               !
950               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
951               DO jk=1,jpkm1             
952                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
953               END DO
954               !
955               ! Update barotropic velocities:
956               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
957                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
958                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
959                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
960                  END IF
961               ENDIF             
962               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
963               !       
964               ! Correct "before" velocities to hold correct bt component:
965               spgv(ji,jj) = 0.e0
966               DO jk=1,jpkm1
967                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
968               END DO
969               !
970               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
971               DO jk=1,jpkm1             
972                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
973               END DO
974               !
975            END DO
976         END DO
977         !
978         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
979            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
980         ENDIF
981         !
982      ENDIF
983      !
984   END SUBROUTINE updatev2d
985
986
987   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
988      !!----------------------------------------------------------------------
989      !!                   *** ROUTINE updateSSH ***
990      !!----------------------------------------------------------------------
991      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
992      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
993      LOGICAL                         , INTENT(in   ) ::   before
994      !!
995      INTEGER :: ji, jj
996      !!----------------------------------------------------------------------
997      !
998      IF( before ) THEN
999         DO jj=j1,j2
1000            DO ji=i1,i2
1001               tabres(ji,jj) = sshn(ji,jj)
1002            END DO
1003         END DO
1004      ELSE
1005         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
1006            DO jj=j1,j2
1007               DO ji=i1,i2
1008                  sshb(ji,jj) =   sshb(ji,jj) &
1009                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
1010               END DO
1011            END DO
1012         ENDIF
1013         !
1014         DO jj=j1,j2
1015            DO ji=i1,i2
1016               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
1017            END DO
1018         END DO
1019         !
1020         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1021            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
1022         ENDIF
1023         !
1024
1025      ENDIF
1026      !
1027   END SUBROUTINE updateSSH
1028
1029
1030   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
1031      !!----------------------------------------------------------------------
1032      !!                      *** ROUTINE updateub2b ***
1033      !!----------------------------------------------------------------------
1034      INTEGER                            , INTENT(in) ::   i1, i2, j1, j2
1035      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
1036      LOGICAL                            , INTENT(in) ::   before
1037      !!
1038      INTEGER :: ji, jj
1039      REAL(wp) :: zrhoy, za1, zcor
1040      !!---------------------------------------------
1041      !
1042      IF (before) THEN
1043         zrhoy = Agrif_Rhoy()
1044         DO jj=j1,j2
1045            DO ji=i1,i2
1046               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1047            END DO
1048         END DO
1049         tabres = zrhoy * tabres
1050      ELSE
1051         !
1052         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1053         !
1054         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1055         DO jj=j1,j2
1056            DO ji=i1,i2
1057               zcor=tabres(ji,jj) - ub2_b(ji,jj)
1058               ! Update time integrated fluxes also in case of multiply nested grids:
1059               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor 
1060               ! Update corrective fluxes:
1061               un_bf(ji,jj)  = un_bf(ji,jj) + zcor
1062               ! Update half step back fluxes:
1063               ub2_b(ji,jj) = tabres(ji,jj)
1064            END DO
1065         END DO
1066      ENDIF
1067      !
1068   END SUBROUTINE updateub2b
1069
1070   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
1071      !!---------------------------------------------
1072      !!          *** ROUTINE reflux_sshu ***
1073      !!---------------------------------------------
1074      INTEGER, INTENT(in) :: i1, i2, j1, j2
1075      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1076      LOGICAL, INTENT(in) :: before
1077      INTEGER, INTENT(in) :: nb, ndir
1078      !!
1079      LOGICAL :: western_side, eastern_side 
1080      INTEGER :: ji, jj
1081      REAL(wp) :: zrhoy, za1, zcor
1082      !!---------------------------------------------
1083      !
1084      IF (before) THEN
1085         zrhoy = Agrif_Rhoy()
1086         DO jj=j1,j2
1087            DO ji=i1,i2
1088               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1089            END DO
1090         END DO
1091         tabres = zrhoy * tabres
1092      ELSE
1093         !
1094         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1095         !
1096         western_side  = (nb == 1).AND.(ndir == 1)
1097         eastern_side  = (nb == 1).AND.(ndir == 2)
1098         !
1099         IF (western_side) THEN
1100            DO jj=j1,j2
1101               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
1102               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
1103               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor
1104            END DO
1105         ENDIF
1106         IF (eastern_side) THEN
1107            DO jj=j1,j2
1108               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
1109               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
1110               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor
1111            END DO
1112         ENDIF
1113         !
1114      ENDIF
1115      !
1116   END SUBROUTINE reflux_sshu
1117
1118   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
1119      !!----------------------------------------------------------------------
1120      !!                      *** ROUTINE updatevb2b ***
1121      !!----------------------------------------------------------------------
1122      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1123      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
1124      LOGICAL                         , INTENT(in   ) ::   before
1125      !!
1126      INTEGER :: ji, jj
1127      REAL(wp) :: zrhox, za1, zcor
1128      !!---------------------------------------------
1129      !
1130      IF( before ) THEN
1131         zrhox = Agrif_Rhox()
1132         DO jj=j1,j2
1133            DO ji=i1,i2
1134               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
1135            END DO
1136         END DO
1137         tabres = zrhox * tabres
1138      ELSE
1139         !
1140         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1141         !
1142         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1143         DO jj=j1,j2
1144            DO ji=i1,i2
1145               zcor=tabres(ji,jj) - vb2_b(ji,jj)
1146               ! Update time integrated fluxes also in case of multiply nested grids:
1147               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor 
1148               ! Update corrective fluxes:
1149               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor
1150               ! Update half step back fluxes:
1151               vb2_b(ji,jj) = tabres(ji,jj)
1152            END DO
1153         END DO
1154      ENDIF
1155      !
1156   END SUBROUTINE updatevb2b
1157
1158   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
1159      !!---------------------------------------------
1160      !!          *** ROUTINE reflux_sshv ***
1161      !!---------------------------------------------
1162      INTEGER, INTENT(in) :: i1, i2, j1, j2
1163      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1164      LOGICAL, INTENT(in) :: before
1165      INTEGER, INTENT(in) :: nb, ndir
1166      !!
1167      LOGICAL :: southern_side, northern_side 
1168      INTEGER :: ji, jj
1169      REAL(wp) :: zrhox, za1, zcor
1170      !!---------------------------------------------
1171      !
1172      IF (before) THEN
1173         zrhox = Agrif_Rhox()
1174         DO jj=j1,j2
1175            DO ji=i1,i2
1176               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
1177            END DO
1178         END DO
1179         tabres = zrhox * tabres
1180      ELSE
1181         !
1182         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1183         !
1184         southern_side = (nb == 2).AND.(ndir == 1)
1185         northern_side = (nb == 2).AND.(ndir == 2)
1186         !
1187         IF (southern_side) THEN
1188            DO ji=i1,i2
1189               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
1190               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
1191               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor
1192            END DO
1193         ENDIF
1194         IF (northern_side) THEN               
1195            DO ji=i1,i2
1196               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
1197               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
1198               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
1199            END DO
1200         ENDIF
1201         !
1202      ENDIF
1203      !
1204   END SUBROUTINE reflux_sshv
1205
1206   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
1207      !
1208      ! ====>>>>>>>>>>    currently not used
1209      !
1210      !!----------------------------------------------------------------------
1211      !!                      *** ROUTINE updateT ***
1212      !!----------------------------------------------------------------------
1213      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
1214      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
1215      LOGICAL                                    , INTENT(in   ) ::   before
1216      !!
1217      INTEGER :: ji,jj,jk
1218      REAL(wp) :: ztemp
1219      !!----------------------------------------------------------------------
1220
1221      IF (before) THEN
1222         DO jk=k1,k2
1223            DO jj=j1,j2
1224               DO ji=i1,i2
1225                  tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
1226                  tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
1227                  tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
1228               END DO
1229            END DO
1230         END DO
1231         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
1232         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
1233         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
1234      ELSE
1235         DO jk=k1,k2
1236            DO jj=j1,j2
1237               DO ji=i1,i2
1238                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
1239                     print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
1240                     print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
1241                     print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
1242                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
1243                     print *,'CORR = ',ztemp-1.
1244                     print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
1245                           tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
1246                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
1247                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
1248                  END IF
1249               END DO
1250            END DO
1251         END DO
1252      ENDIF
1253      !
1254   END SUBROUTINE update_scales
1255
1256
1257   SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
1258      !!----------------------------------------------------------------------
1259      !!                      *** ROUTINE updateen ***
1260      !!----------------------------------------------------------------------
1261      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1262      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1263      LOGICAL                               , INTENT(in   ) ::   before
1264      !!----------------------------------------------------------------------
1265      !
1266      IF( before ) THEN
1267         ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
1268      ELSE
1269         en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1270      ENDIF
1271      !
1272   END SUBROUTINE updateEN
1273
1274
1275   SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
1276      !!----------------------------------------------------------------------
1277      !!                      *** ROUTINE updateavt ***
1278      !!----------------------------------------------------------------------
1279      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1280      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1281      LOGICAL                               , INTENT(in   ) ::   before
1282      !!----------------------------------------------------------------------
1283      !
1284      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
1285      ELSE                ;   avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1286      ENDIF
1287      !
1288   END SUBROUTINE updateAVT
1289
1290
1291   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
1292      !!---------------------------------------------
1293      !!           *** ROUTINE updateavm ***
1294      !!----------------------------------------------------------------------
1295      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1296      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1297      LOGICAL                               , INTENT(in   ) ::   before
1298      !!----------------------------------------------------------------------
1299      !
1300      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1301      ELSE                ;   avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1302      ENDIF
1303      !
1304   END SUBROUTINE updateAVM
1305
1306   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
1307      !!---------------------------------------------
1308      !!           *** ROUTINE updatee3t ***
1309      !!---------------------------------------------
1310      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1311      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1312      LOGICAL, INTENT(in) :: before
1313      !
1314      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1315      INTEGER :: ji,jj,jk
1316      REAL(wp) :: zcoef
1317      !!---------------------------------------------
1318      !
1319      IF (.NOT.before) THEN
1320         !
1321         ALLOCATE(ptab(i1:i2,j1:j2,1:jpk)) 
1322         !
1323         ! Update e3t from ssh (z* case only)
1324         DO jk = 1, jpkm1
1325            DO jj=j1,j2
1326               DO ji=i1,i2
1327                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
1328                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1329               END DO
1330            END DO
1331         END DO
1332         !
1333         ! 1) Updates at BEFORE time step:
1334         ! -------------------------------
1335         !
1336         ! Save "old" scale factor (prior update) for subsequent asselin correction
1337         ! of prognostic variables
1338         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
1339
1340         ! One should also save e3t_b, but lacking of workspace...
1341!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1)
1342
1343         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
1344            DO jk = 1, jpkm1
1345               DO jj=j1,j2
1346                  DO ji=i1,i2
1347                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
1348                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
1349                  END DO
1350               END DO
1351            END DO
1352            !
1353            e3w_b  (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
1354            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
1355            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
1356            !
1357            DO jk = 2, jpk
1358               DO jj = j1,j2
1359                  DO ji = i1,i2           
1360                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1361                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
1362                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
1363                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
1364                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1365                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
1366                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
1367                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
1368                  END DO
1369               END DO
1370            END DO
1371            !
1372         ENDIF       
1373         !
1374         ! 2) Updates at NOW time step:
1375         ! ----------------------------
1376         !
1377         ! Update vertical scale factor at T-points:
1378         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
1379         !
1380         ! Update total depth:
1381         ht_n(i1:i2,j1:j2) = 0._wp
1382         DO jk = 1, jpkm1
1383            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk)
1384         END DO
1385         !
1386         ! Update vertical scale factor at W-points and depths:
1387         e3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
1388         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
1389         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
1390         gde3w_n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
1391         !
1392         DO jk = 2, jpk
1393            DO jj = j1,j2
1394               DO ji = i1,i2           
1395               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1396               e3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   &
1397               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1398               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
1399               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
1400                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
1401               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
1402               END DO
1403            END DO
1404         END DO
1405         !
1406         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1407            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
1408            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
1409            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
1410            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
1411         ENDIF
1412         !
1413         DEALLOCATE(ptab)
1414      ENDIF
1415      !
1416   END SUBROUTINE updatee3t
1417
1418#else
1419   !!----------------------------------------------------------------------
1420   !!   Empty module                                          no AGRIF zoom
1421   !!----------------------------------------------------------------------
1422CONTAINS
1423   SUBROUTINE agrif_oce_update_empty
1424      WRITE(*,*)  'agrif_oce_update : You should not have seen this print! error?'
1425   END SUBROUTINE agrif_oce_update_empty
1426#endif
1427
1428   !!======================================================================
1429END MODULE agrif_oce_update
1430
Note: See TracBrowser for help on using the repository browser.