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/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/NST – NEMO

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/NST/agrif_oce_update.F90 @ 12928

Last change on this file since 12928 was 12928, checked in by smueller, 4 years ago

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

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