source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_update.F90 @ 13351

Last change on this file since 13351 was 13351, checked in by jchanut, 7 months ago

#2222, fixes mistakes when removing key_vertical - rename flag in namelist to prevent further mistakes - vertical refinement works fine in VORTEX

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