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

Last change on this file was 15265, checked in by jchanut, 7 days ago

Improve ghost cell initialization with AGRIF + minor changes such as missing _wp, tests namelists updates, etc… can be assigned to #2638

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