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/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_oce_update.F90 @ 14064

Last change on this file since 14064 was 14064, checked in by ayoung, 22 months ago

Merging ticket #2506 into trunk.

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