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_opa_update.F90 in branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC – NEMO

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

Last change on this file since 8741 was 8741, checked in by jchanut, 5 years ago

AGRIF + vvl Main changes - #1965

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