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_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 8010

Last change on this file since 8010 was 8010, checked in by jchanut, 7 years ago

AGRIF vvl add on

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