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

Last change on this file since 14086 was 14086, checked in by cetlod, 22 months ago

Adding AGRIF branches into the trunk

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