New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_oce_update.F90 in NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/NST – NEMO

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/NST/agrif_oce_update.F90 @ 10419

Last change on this file since 10419 was 10068, checked in by nicolasmartin, 6 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

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