source: branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 8965

Last change on this file since 8965 was 8965, checked in by jchanut, 3 years ago

Add zoom capability in default Gyre, remove bcs velocity update (non reproducible), add option to set clamped obcs - #1965

  • Property svn:keywords set to Id
File size: 41.3 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_opa_update
6#if defined key_agrif 
7   USE par_oce
8   USE oce
9   USE dom_oce
10   USE agrif_oce
11   USE in_out_manager  ! I/O manager
12   USE lib_mpp
13   USE wrk_nemo 
14   USE zdf_oce        ! vertical physics: ocean variables
15   USE domvvl         ! Need interpolation routines
16
17   IMPLICIT NONE
18   PRIVATE
19
20   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales, Agrif_Update_vvl, Agrif_Update_ssh
21
22# if defined key_zdftke
23   PUBLIC Agrif_Update_Tke
24# endif
25   !!----------------------------------------------------------------------
26   !! NEMO/NST 3.6 , NEMO Consortium (2010)
27   !! $Id$
28   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
29   !!----------------------------------------------------------------------
30CONTAINS
31
32   SUBROUTINE Agrif_Update_Tra( )
33      !!---------------------------------------------
34      !!   *** ROUTINE Agrif_Update_Tra ***
35      !!---------------------------------------------
36      !
37      IF (Agrif_Root()) RETURN
38      !
39#if defined TWO_WAY 
40      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed()
41
42      Agrif_UseSpecialValueInUpdate = .TRUE.
43      Agrif_SpecialValueFineGrid = 0.
44      !
45# if ! defined DECAL_FEEDBACK
46      CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
47! near boundary update:
48!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
49# else
50      CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
51! near boundary update:
52!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
53# endif
54      !
55      Agrif_UseSpecialValueInUpdate = .FALSE.
56      !
57#endif
58      !
59   END SUBROUTINE Agrif_Update_Tra
60
61   SUBROUTINE Agrif_Update_Dyn( )
62      !!---------------------------------------------
63      !!   *** ROUTINE Agrif_Update_Dyn ***
64      !!---------------------------------------------
65      !
66      IF (Agrif_Root()) RETURN
67      !
68#if defined TWO_WAY
69      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed()
70
71      Agrif_UseSpecialValueInUpdate = .FALSE.
72      Agrif_SpecialValueFineGrid = 0.
73      !     
74# if ! defined DECAL_FEEDBACK
75      CALL Agrif_Update_Variable(un_update_id,procname = updateU)
76      CALL Agrif_Update_Variable(vn_update_id,procname = updateV)
77! near boundary update:
78!      CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
79!      CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)
80# else
81      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
82      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
83! near boundary update:
84!      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
85!      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
86# endif
87
88# if ! defined DECAL_FEEDBACK
89      CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
90      CALL Agrif_Update_Variable(e2v_id,procname = updateV2d) 
91# else
92      CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
93      CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d) 
94# endif
95      !
96# if ! defined DECAL_FEEDBACK
97      ! Account for updated thicknesses at boundary edges
98      IF (.NOT.ln_linssh) THEN
99! For the time being calls below do not ensure reproducible results
100!         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)
101!         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)
102      ENDIF
103# endif
104      !
105      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
106         ! Update time integrated transports
107#  if ! defined DECAL_FEEDBACK
108         CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
109         CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
110#  else
111         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
112         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
113#  endif
114      END IF
115#endif
116      !
117   END SUBROUTINE Agrif_Update_Dyn
118
119   SUBROUTINE Agrif_Update_ssh( )
120      !!---------------------------------------------
121      !!   *** ROUTINE Agrif_Update_ssh ***
122      !!---------------------------------------------
123      !
124      IF (Agrif_Root()) RETURN
125      !
126#if defined TWO_WAY
127      !
128      Agrif_UseSpecialValueInUpdate = .TRUE.
129      Agrif_SpecialValueFineGrid = 0.
130# if ! defined DECAL_FEEDBACK
131      CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
132# else
133      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
134# endif
135      !
136      Agrif_UseSpecialValueInUpdate = .FALSE.
137      !
138#  if defined DECAL_FEEDBACK && defined VOL_REFLUX
139      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
140         ! Refluxing on ssh:
141         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)
142         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv)
143      END IF
144#  endif
145      !
146#endif
147      !
148   END SUBROUTINE Agrif_Update_ssh
149
150# if defined key_zdftke
151
152   SUBROUTINE Agrif_Update_Tke( kt )
153      !!---------------------------------------------
154      !!   *** ROUTINE Agrif_Update_Tke ***
155      !!---------------------------------------------
156      !!
157      INTEGER, INTENT(in) :: kt 
158      !
159      IF (Agrif_Root()) RETURN
160      !       
161      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
162#  if defined TWO_WAY
163
164      Agrif_UseSpecialValueInUpdate = .TRUE.
165      Agrif_SpecialValueFineGrid = 0.
166
167      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  )
168      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
169      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
170
171      Agrif_UseSpecialValueInUpdate = .FALSE.
172
173#  endif
174     
175   END SUBROUTINE Agrif_Update_Tke
176   
177# endif /* key_zdftke */
178
179   SUBROUTINE Agrif_Update_vvl( )
180      !!---------------------------------------------
181      !!   *** ROUTINE Agrif_Update_vvl ***
182      !!---------------------------------------------
183      !
184      IF (Agrif_Root()) RETURN
185      !
186#if defined TWO_WAY 
187      !
188      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step()
189      !
190      Agrif_UseSpecialValueInUpdate = .TRUE.
191      Agrif_SpecialValueFineGrid = 0.
192      !
193      ! No interface separation here, update vertical grid at T points
194      ! everywhere over the overlapping regions (one account for refluxing in that case):
195      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t) 
196      !
197      Agrif_UseSpecialValueInUpdate = .FALSE.
198      !
199      CALL Agrif_ChildGrid_To_ParentGrid()
200      CALL dom_vvl_update_UVF
201      CALL Agrif_ParentGrid_To_ChildGrid()
202      !
203#endif
204      !
205   END SUBROUTINE Agrif_Update_vvl
206
207   SUBROUTINE dom_vvl_update_UVF
208      !!---------------------------------------------
209      !!       *** ROUTINE dom_vvl_update_UVF ***
210      !!---------------------------------------------
211      !!
212      INTEGER :: jk
213      REAL(wp):: zcoef
214      !!---------------------------------------------
215
216      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', &
217                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step()
218
219      ! Save "old" scale factor (prior update) for subsequent asselin correction
220      ! of prognostic variables
221      ! -----------------------
222      !
223      e3u_a(:,:,:) = e3u_n(:,:,:)
224      e3v_a(:,:,:) = e3v_n(:,:,:)
225!      ua(:,:,:) = e3u_b(:,:,:)
226!      va(:,:,:) = e3v_b(:,:,:)
227      hu_a(:,:) = hu_n(:,:)
228      hv_a(:,:) = hv_n(:,:)
229
230      ! 1) NOW fields
231      !--------------
232     
233         ! Vertical scale factor interpolations
234         ! ------------------------------------
235      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' )
236      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' )
237      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' )
238
239      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
240      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
241
242         ! Update total depths:
243         ! --------------------
244      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points
245      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points
246      DO jk = 1, jpkm1
247         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
248         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
249      END DO
250      !                                        ! Inverse of the local depth
251      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
252      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
253
254
255      ! 2) BEFORE fields:
256      !------------------
257!      IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
258!         & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
259!         & .AND.(.NOT.ln_bt_fw)))) THEN
260      IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
261         !
262         ! Vertical scale factor interpolations
263         ! ------------------------------------
264         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  )
265         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  )
266
267         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
268         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
269
270         ! Update total depths:
271         ! --------------------
272         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points
273         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points
274         DO jk = 1, jpkm1
275            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
276            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
277         END DO
278         !                                     ! Inverse of the local depth
279         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
280         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
281      ENDIF
282      !
283   END SUBROUTINE dom_vvl_update_UVF
284
285   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
286      !!---------------------------------------------
287      !!           *** ROUTINE updateT ***
288      !!---------------------------------------------
289      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
290      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
291      LOGICAL, INTENT(in) :: before
292      !!
293      INTEGER :: ji,jj,jk,jn
294      REAL(wp) :: ztb, ztnu, ztno
295      !!---------------------------------------------
296      !
297      IF (before) THEN
298         DO jn = n1,n2
299            DO jk=k1,k2
300               DO jj=j1,j2
301                  DO ji=i1,i2
302!> jc tmp
303                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
304!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
305!< jc tmp
306                  END DO
307               END DO
308            END DO
309         END DO
310      ELSE
311!> jc tmp
312         DO jn = n1,n2
313            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
314                                         & * tmask(i1:i2,j1:j2,k1:k2)
315         ENDDO
316!< jc tmp
317         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
318            ! Add asselin part
319            DO jn = n1,n2
320               DO jk=k1,k2
321                  DO jj=j1,j2
322                     DO ji=i1,i2
323                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
324                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
325                           ztnu = tabres(ji,jj,jk,jn)
326                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
327                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
328                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
329                        ENDIF
330                     ENDDO
331                  ENDDO
332               ENDDO
333            ENDDO
334         ENDIF
335         DO jn = n1,n2
336            DO jk=k1,k2
337               DO jj=j1,j2
338                  DO ji=i1,i2
339                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
340                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
341                     END IF
342                  END DO
343               END DO
344            END DO
345         END DO
346         !
347         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
348            tsb(i1:i2,j1:j2,k1:k2,n1:n2)  = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
349         ENDIF
350         !
351      ENDIF
352      !
353   END SUBROUTINE updateTS
354
355
356   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before )
357      !!---------------------------------------------
358      !!           *** ROUTINE updateu ***
359      !!---------------------------------------------
360      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
361      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
362      LOGICAL                               , INTENT(in   ) :: before
363      !
364      INTEGER  :: ji, jj, jk
365      REAL(wp) :: zrhoy, zub, zunu, zuno
366      !!---------------------------------------------
367      !
368      IF( before ) THEN
369         zrhoy = Agrif_Rhoy()
370         DO jk = k1, k2
371            tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
372         END DO
373      ELSE
374         DO jk=k1,k2
375            DO jj=j1,j2
376               DO ji=i1,i2
377                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) 
378                  !
379                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
380                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used
381                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
382                     zunu = tabres(ji,jj,jk)
383                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &     
384                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
385                  ENDIF
386                  !
387                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
388               END DO
389            END DO
390         END DO
391         !
392         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
393            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
394         ENDIF
395         !
396      ENDIF
397      !
398   END SUBROUTINE updateu
399
400   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir )
401      !!---------------------------------------------
402      !!           *** ROUTINE correct_u_bdy ***
403      !!---------------------------------------------
404      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
405      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
406      LOGICAL                               , INTENT(in   ) :: before
407      INTEGER                               , INTENT(in)    :: nb, ndir
408      !!
409      LOGICAL :: western_side, eastern_side 
410      !
411      INTEGER  :: jj, jk
412      REAL(wp) :: zcor
413      !!---------------------------------------------
414      !
415      IF( .NOT.before ) THEN
416         !
417         western_side  = (nb == 1).AND.(ndir == 1)
418         eastern_side  = (nb == 1).AND.(ndir == 2)
419         !
420         IF (western_side) THEN
421            DO jj=j1,j2
422               zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
423               un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
424               DO jk=1,jpkm1
425                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)
426               END DO
427            END DO
428         ENDIF
429         !
430         IF (eastern_side) THEN
431            DO jj=j1,j2
432               zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
433               un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
434               DO jk=1,jpkm1
435                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)
436               END DO
437            END DO
438         ENDIF
439         !
440      ENDIF
441      !
442   END SUBROUTINE correct_u_bdy
443
444
445   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before)
446      !!---------------------------------------------
447      !!           *** ROUTINE updatev ***
448      !!---------------------------------------------
449      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
450      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
451      LOGICAL                               , INTENT(in   ) :: before
452      !
453      INTEGER  :: ji, jj, jk
454      REAL(wp) :: zrhox, zvb, zvnu, zvno
455      !!---------------------------------------------     
456      !
457      IF (before) THEN
458         zrhox = Agrif_Rhox()
459         DO jk=k1,k2
460            DO jj=j1,j2
461               DO ji=i1,i2
462                  tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
463               END DO
464            END DO
465         END DO
466      ELSE
467         DO jk=k1,k2
468            DO jj=j1,j2
469               DO ji=i1,i2
470                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj)
471                  !
472                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
473                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
474                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
475                     zvnu = tabres(ji,jj,jk)
476                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &     
477                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
478                  ENDIF
479                  !
480                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
481               END DO
482            END DO
483         END DO
484         !
485         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
486            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
487         ENDIF
488         !
489      ENDIF
490      !
491   END SUBROUTINE updatev
492
493   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir )
494      !!---------------------------------------------
495      !!           *** ROUTINE correct_u_bdy ***
496      !!---------------------------------------------
497      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
498      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
499      LOGICAL                               , INTENT(in   ) :: before
500      INTEGER                               , INTENT(in)    :: nb, ndir
501      !!
502      LOGICAL :: southern_side, northern_side 
503      !
504      INTEGER  :: ji, jk
505      REAL(wp) :: zcor
506      !!---------------------------------------------
507      !
508      IF( .NOT.before ) THEN
509         !
510         southern_side = (nb == 2).AND.(ndir == 1)
511         northern_side = (nb == 2).AND.(ndir == 2)
512         !
513         IF (southern_side) THEN
514            DO ji=i1,i2
515               zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
516               vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
517               DO jk=1,jpkm1
518                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)
519               END DO
520            END DO
521         ENDIF
522         !
523         IF (northern_side) THEN
524            DO ji=i1,i2
525               zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
526               vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
527               DO jk=1,jpkm1
528                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)
529               END DO
530            END DO
531         ENDIF
532         !
533      ENDIF
534      !
535   END SUBROUTINE correct_v_bdy
536
537
538   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
539      !!---------------------------------------------
540      !!          *** ROUTINE updateu2d ***
541      !!---------------------------------------------
542      INTEGER, INTENT(in) :: i1, i2, j1, j2
543      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
544      LOGICAL, INTENT(in) :: before
545      !!
546      INTEGER  :: ji, jj, jk
547      REAL(wp) :: zrhoy
548      REAL(wp) :: zcorr
549      !!---------------------------------------------
550      !
551      IF (before) THEN
552         zrhoy = Agrif_Rhoy()
553         DO jj=j1,j2
554            DO ji=i1,i2
555               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
556            END DO
557         END DO
558      ELSE
559         DO jj=j1,j2
560            DO ji=i1,i2
561               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj) 
562               !   
563               ! Update "now" 3d velocities:
564               spgu(ji,jj) = 0._wp
565               DO jk=1,jpkm1
566                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
567               END DO
568               !
569               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
570               DO jk=1,jpkm1             
571                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
572               END DO
573               !
574               ! Update barotropic velocities:
575               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
576                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
577                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
578                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
579                  END IF
580               ENDIF   
581               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
582               !       
583               ! Correct "before" velocities to hold correct bt component:
584               spgu(ji,jj) = 0.e0
585               DO jk=1,jpkm1
586                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
587               END DO
588               !
589               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
590               DO jk=1,jpkm1             
591                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
592               END DO
593               !
594            END DO
595         END DO
596         !
597         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
598            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
599         ENDIF
600      ENDIF
601      !
602   END SUBROUTINE updateu2d
603
604
605   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
606      !!---------------------------------------------
607      !!          *** ROUTINE updatev2d ***
608      !!---------------------------------------------
609      INTEGER, INTENT(in) :: i1, i2, j1, j2
610      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
611      LOGICAL, INTENT(in) :: before
612      !!
613      INTEGER  :: ji, jj, jk
614      REAL(wp) :: zrhox
615      REAL(wp) :: zcorr
616      !!---------------------------------------------
617      !
618      IF (before) THEN
619         zrhox = Agrif_Rhox()
620         DO jj=j1,j2
621            DO ji=i1,i2
622               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
623            END DO
624         END DO
625      ELSE
626         DO jj=j1,j2
627            DO ji=i1,i2
628               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj) 
629               !   
630               ! Update "now" 3d velocities:
631               spgv(ji,jj) = 0.e0
632               DO jk=1,jpkm1
633                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
634               END DO
635               !
636               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
637               DO jk=1,jpkm1             
638                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
639               END DO
640               !
641               ! Update barotropic velocities:
642               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
643                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
644                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
645                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
646                  END IF
647               ENDIF             
648               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
649               !       
650               ! Correct "before" velocities to hold correct bt component:
651               spgv(ji,jj) = 0.e0
652               DO jk=1,jpkm1
653                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
654               END DO
655               !
656               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
657               DO jk=1,jpkm1             
658                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
659               END DO
660               !
661            END DO
662         END DO
663         !
664         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
665            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
666         ENDIF
667         !
668      ENDIF
669      !
670   END SUBROUTINE updatev2d
671
672
673   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
674      !!---------------------------------------------
675      !!          *** ROUTINE updateSSH ***
676      !!---------------------------------------------
677      INTEGER, INTENT(in) :: i1, i2, j1, j2
678      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
679      LOGICAL, INTENT(in) :: before
680      !!
681      INTEGER :: ji, jj
682      !!---------------------------------------------
683      !
684      IF (before) THEN
685         DO jj=j1,j2
686            DO ji=i1,i2
687               tabres(ji,jj) = sshn(ji,jj)
688            END DO
689         END DO
690      ELSE
691         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
692            DO jj=j1,j2
693               DO ji=i1,i2
694                  sshb(ji,jj) =   sshb(ji,jj) &
695                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
696               END DO
697            END DO
698         ENDIF
699         !
700         DO jj=j1,j2
701            DO ji=i1,i2
702               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
703            END DO
704         END DO
705         !
706         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
707            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
708         ENDIF
709         !
710
711      ENDIF
712      !
713   END SUBROUTINE updateSSH
714
715
716   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
717      !!---------------------------------------------
718      !!          *** ROUTINE updateub2b ***
719      !!---------------------------------------------
720      INTEGER, INTENT(in) :: i1, i2, j1, j2
721      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
722      LOGICAL, INTENT(in) :: before
723      !!
724      INTEGER :: ji, jj
725      REAL(wp) :: zrhoy, za1, zcor
726      !!---------------------------------------------
727      !
728      IF (before) THEN
729         zrhoy = Agrif_Rhoy()
730         DO jj=j1,j2
731            DO ji=i1,i2
732               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
733            END DO
734         END DO
735         tabres = zrhoy * tabres
736      ELSE
737         !
738         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
739         !
740         za1 = 1._wp / REAL(Agrif_rhot(), wp)
741         DO jj=j1,j2
742            DO ji=i1,i2
743               zcor=tabres(ji,jj) - ub2_b(ji,jj)
744               ! Update time integrated fluxes also in case of multiply nested grids:
745               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor 
746               ! Update corrective fluxes:
747               un_bf(ji,jj)  = un_bf(ji,jj) + zcor
748               ! Update half step back fluxes:
749               ub2_b(ji,jj) = tabres(ji,jj)
750            END DO
751         END DO
752      ENDIF
753      !
754   END SUBROUTINE updateub2b
755
756   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
757      !!---------------------------------------------
758      !!          *** ROUTINE reflux_sshu ***
759      !!---------------------------------------------
760      INTEGER, INTENT(in) :: i1, i2, j1, j2
761      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
762      LOGICAL, INTENT(in) :: before
763      INTEGER, INTENT(in) :: nb, ndir
764      !!
765      LOGICAL :: western_side, eastern_side 
766      INTEGER :: ji, jj
767      REAL(wp) :: zrhoy, za1, zcor
768      !!---------------------------------------------
769      !
770      IF (before) THEN
771         zrhoy = Agrif_Rhoy()
772         DO jj=j1,j2
773            DO ji=i1,i2
774               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
775            END DO
776         END DO
777         tabres = zrhoy * tabres
778      ELSE
779         !
780         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
781         !
782         western_side  = (nb == 1).AND.(ndir == 1)
783         eastern_side  = (nb == 1).AND.(ndir == 2)
784         !
785         IF (western_side) THEN
786            DO jj=j1,j2
787               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
788               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
789               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor
790            END DO
791         ENDIF
792         IF (eastern_side) THEN
793            DO jj=j1,j2
794               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
795               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
796               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor
797            END DO
798         ENDIF
799         !
800      ENDIF
801      !
802   END SUBROUTINE reflux_sshu
803
804   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
805      !!---------------------------------------------
806      !!          *** ROUTINE updatevb2b ***
807      !!---------------------------------------------
808      INTEGER, INTENT(in) :: i1, i2, j1, j2
809      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
810      LOGICAL, INTENT(in) :: before
811      !!
812      INTEGER :: ji, jj
813      REAL(wp) :: zrhox, za1, zcor
814      !!---------------------------------------------
815      !
816      IF (before) THEN
817         zrhox = Agrif_Rhox()
818         DO jj=j1,j2
819            DO ji=i1,i2
820               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
821            END DO
822         END DO
823         tabres = zrhox * tabres
824      ELSE
825         !
826         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
827         !
828         za1 = 1._wp / REAL(Agrif_rhot(), wp)
829         DO jj=j1,j2
830            DO ji=i1,i2
831               zcor=tabres(ji,jj) - vb2_b(ji,jj)
832               ! Update time integrated fluxes also in case of multiply nested grids:
833               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor 
834               ! Update corrective fluxes:
835               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor
836               ! Update half step back fluxes:
837               vb2_b(ji,jj) = tabres(ji,jj)
838            END DO
839         END DO
840      ENDIF
841      !
842   END SUBROUTINE updatevb2b
843
844   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
845      !!---------------------------------------------
846      !!          *** ROUTINE reflux_sshv ***
847      !!---------------------------------------------
848      INTEGER, INTENT(in) :: i1, i2, j1, j2
849      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
850      LOGICAL, INTENT(in) :: before
851      INTEGER, INTENT(in) :: nb, ndir
852      !!
853      LOGICAL :: southern_side, northern_side 
854      INTEGER :: ji, jj
855      REAL(wp) :: zrhox, za1, zcor
856      !!---------------------------------------------
857      !
858      IF (before) THEN
859         zrhox = Agrif_Rhox()
860         DO jj=j1,j2
861            DO ji=i1,i2
862               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
863            END DO
864         END DO
865         tabres = zrhox * tabres
866      ELSE
867         !
868         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
869         !
870         southern_side = (nb == 2).AND.(ndir == 1)
871         northern_side = (nb == 2).AND.(ndir == 2)
872         !
873         IF (southern_side) THEN
874            DO ji=i1,i2
875               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
876               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
877               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor
878            END DO
879         ENDIF
880         IF (northern_side) THEN               
881            DO ji=i1,i2
882               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
883               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
884               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
885            END DO
886         ENDIF
887         !
888      ENDIF
889      !
890   END SUBROUTINE reflux_sshv
891
892   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
893      ! currently not used
894      !!---------------------------------------------
895      !!           *** ROUTINE updateT ***
896      !!---------------------------------------------
897      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
898      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
899      LOGICAL, iNTENT(in) :: before
900      !
901      INTEGER :: ji,jj,jk
902      REAL(wp) :: ztemp
903      !!---------------------------------------------
904
905      IF (before) THEN
906         DO jk=k1,k2
907            DO jj=j1,j2
908               DO ji=i1,i2
909                  tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
910                  tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
911                  tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
912               END DO
913            END DO
914         END DO
915         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
916         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
917         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
918      ELSE
919         DO jk=k1,k2
920            DO jj=j1,j2
921               DO ji=i1,i2
922                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
923                     print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
924                     print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
925                     print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
926                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
927                     print *,'CORR = ',ztemp-1.
928                     print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
929                           tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
930                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
931                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
932                  END IF
933               END DO
934            END DO
935         END DO
936      ENDIF
937      !
938   END SUBROUTINE update_scales
939
940# if defined key_zdftke
941
942   SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
943      !!---------------------------------------------
944      !!           *** ROUTINE updateen ***
945      !!---------------------------------------------
946      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
947      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
948      LOGICAL, INTENT(in) :: before
949      !!---------------------------------------------
950      !
951      IF (before) THEN
952         ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
953      ELSE
954         en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
955      ENDIF
956      !
957   END SUBROUTINE updateEN
958
959
960   SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
961      !!---------------------------------------------
962      !!           *** ROUTINE updateavt ***
963      !!---------------------------------------------
964      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
965      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
966      LOGICAL, INTENT(in) :: before
967      !!---------------------------------------------
968      !
969      IF (before) THEN
970         ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
971      ELSE
972         avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
973      ENDIF
974      !
975   END SUBROUTINE updateAVT
976
977
978   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
979      !!---------------------------------------------
980      !!           *** ROUTINE updateavm ***
981      !!---------------------------------------------
982      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
983      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
984      LOGICAL, INTENT(in) :: before
985      !!---------------------------------------------
986      !
987      IF (before) THEN
988         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
989      ELSE
990         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
991      ENDIF
992      !
993   END SUBROUTINE updateAVM
994
995# endif /* key_zdftke */ 
996
997   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
998      !!---------------------------------------------
999      !!           *** ROUTINE updatee3t ***
1000      !!---------------------------------------------
1001      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1002      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1003      LOGICAL, INTENT(in) :: before
1004      !
1005      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1006      INTEGER :: ji,jj,jk
1007      REAL(wp) :: zcoef
1008      !!---------------------------------------------
1009      !
1010      IF (.NOT.before) THEN
1011         !
1012         ALLOCATE(ptab(i1:i2,j1:j2,1:jpk)) 
1013         !
1014         ! Update e3t from ssh (z* case only)
1015         DO jk = 1, jpkm1
1016            DO jj=j1,j2
1017               DO ji=i1,i2
1018                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
1019                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1020               END DO
1021            END DO
1022         END DO
1023         !
1024         ! 1) Updates at BEFORE time step:
1025         ! -------------------------------
1026         !
1027         ! Save "old" scale factor (prior update) for subsequent asselin correction
1028         ! of prognostic variables
1029         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
1030
1031         ! One should also save e3t_b, but lacking of workspace...
1032!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1)
1033
1034         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
1035            DO jk = 1, jpkm1
1036               DO jj=j1,j2
1037                  DO ji=i1,i2
1038                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
1039                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
1040                  END DO
1041               END DO
1042            END DO
1043            !
1044            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)
1045            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
1046            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
1047            !
1048            DO jk = 2, jpk
1049               DO jj = j1,j2
1050                  DO ji = i1,i2           
1051                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1052                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
1053                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
1054                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
1055                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1056                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
1057                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
1058                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
1059                  END DO
1060               END DO
1061            END DO
1062            !
1063         ENDIF       
1064         !
1065         ! 2) Updates at NOW time step:
1066         ! ----------------------------
1067         !
1068         ! Update vertical scale factor at T-points:
1069         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
1070         !
1071         ! Update total depth:
1072         ht_n(i1:i2,j1:j2) = 0._wp
1073         DO jk = 1, jpkm1
1074            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)
1075         END DO
1076         !
1077         ! Update vertical scale factor at W-points and depths:
1078         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)
1079         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
1080         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
1081         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
1082         !
1083         DO jk = 2, jpk
1084            DO jj = j1,j2
1085               DO ji = i1,i2           
1086               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1087               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) )   &
1088               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1089               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
1090               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
1091                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
1092               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
1093               END DO
1094            END DO
1095         END DO
1096         !
1097         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1098            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
1099            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
1100            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
1101            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
1102         ENDIF
1103         !
1104         DEALLOCATE(ptab)
1105      ENDIF
1106      !
1107   END SUBROUTINE updatee3t
1108
1109#else
1110CONTAINS
1111   SUBROUTINE agrif_opa_update_empty
1112      !!---------------------------------------------
1113      !!   *** ROUTINE agrif_opa_update_empty ***
1114      !!---------------------------------------------
1115      WRITE(*,*)  'agrif_opa_update : You should not have seen this print! error?'
1116   END SUBROUTINE agrif_opa_update_empty
1117#endif
1118END MODULE agrif_opa_update
1119
Note: See TracBrowser for help on using the repository browser.