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

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

AGRIF: Remove update frequency parameter from namelist - #1965

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