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

Last change on this file since 11574 was 11574, checked in by jchanut, 23 months ago

#2222, import changes from dev_r10973_AGRIF-01_jchanut_small_jpi_jpj (i.e. #2199)

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