source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90 @ 11590

Last change on this file since 11590 was 11590, checked in by jchanut, 16 months ago

#2222: 1) create remapping module (vremap) and integration of D. Engwirda piecewise polynomial recontruction package (PPR_LIB cpp key). 2) Various bug corrections with key_vertical activated.

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