source: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 8998

Last change on this file since 8998 was 8998, checked in by timgraham, 3 years ago

First commit of Jerome's modified versions of agrif_opa routines

  • Property svn:keywords set to Id
File size: 53.0 KB
Line 
1#define TWO_WAY        /* TWO WAY NESTING */
2#define 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#if defined key_vertical
285
286   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
287      !!---------------------------------------------
288      !!           *** ROUTINE updateT ***
289      !!---------------------------------------------
290      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
291      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
292      LOGICAL, INTENT(in) :: before
293      !!
294      INTEGER :: ji,jj,jk,jn
295      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
296      REAL(wp) :: h_in(k1:k2)
297      REAL(wp) :: h_out(1:jpk)
298      INTEGER  :: N_in, N_out
299      REAL(wp) :: h_diff
300      REAL(wp) :: zrho_xy
301      REAL(wp) :: tabin(k1:k2,n1:n2)
302      !!---------------------------------------------
303      !
304      IF (before) THEN
305         AGRIF_SpecialValue = -999._wp
306         zrho_xy = Agrif_rhox() * Agrif_rhoy() 
307         DO jn = n1,n2-1
308            DO jk=k1,k2
309               DO jj=j1,j2
310                  DO ji=i1,i2
311                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) &
312                                           * tmask(ji,jj,jk) * zrho_xy + (tmask(ji,jj,jk)-1)*999._wp
313                  END DO
314               END DO
315            END DO
316         END DO
317         DO jk=k1,k2
318            DO jj=j1,j2
319               DO ji=i1,i2
320                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy  &
321                                           + (tmask(ji,jj,jk)-1)*999._wp
322               END DO
323            END DO
324         END DO
325      ELSE
326         tabres_child(:,:,:,:) = 0.
327         AGRIF_SpecialValue = 0._wp
328         DO jj=j1,j2
329            DO ji=i1,i2
330               N_in = 0
331               DO jk=k1,k2 !k2 = jpk of child grid
332                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT
333                  N_in = N_in + 1
334                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
335                  h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj)
336               ENDDO
337               N_out = 0
338               DO jk=1,jpk ! jpk of parent grid
339                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
340                  N_out = N_out + 1
341                  h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
342               ENDDO
343               IF (N_in > 0) THEN !Remove this?
344                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
345                  IF (h_diff < -1.e-4) THEN
346                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
347                     print *,h_in(1:N_in)
348                     print *,h_out(1:N_out)
349                     STOP
350                  ENDIF
351                  DO jn=n1,n2-1
352                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
353                  ENDDO
354               ENDIF
355            ENDDO
356         ENDDO
357
358         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
359            ! Add asselin part
360            DO jn = n1,n2-1
361               DO jk=1,jpk
362                  DO jj=j1,j2
363                     DO ji=i1,i2
364                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
365                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
366                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
367                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
368                        ENDIF
369                     ENDDO
370                  ENDDO
371               ENDDO
372            ENDDO
373         ENDIF
374         DO jn = n1,n2-1
375            DO jk=1,jpk
376               DO jj=j1,j2
377                  DO ji=i1,i2
378                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
379                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
380                     END IF
381                  END DO
382               END DO
383            END DO
384         END DO
385      ENDIF
386      !
387   END SUBROUTINE updateTS
388
389# else
390
391   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
392      !!---------------------------------------------
393      !!           *** ROUTINE updateT ***
394      !!---------------------------------------------
395      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
396      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
397      LOGICAL, INTENT(in) :: before
398      !!
399      INTEGER :: ji,jj,jk,jn
400      REAL(wp) :: ztb, ztnu, ztno
401      !!---------------------------------------------
402      !
403      IF (before) THEN
404         DO jn = 1,jpts
405            DO jk=k1,k2
406               DO jj=j1,j2
407                  DO ji=i1,i2
408!> jc tmp
409                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
410!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
411!< jc tmp
412                  END DO
413               END DO
414            END DO
415         END DO
416      ELSE
417!> jc tmp
418         DO jn = 1,jpts
419            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
420                                         & * tmask(i1:i2,j1:j2,k1:k2)
421         ENDDO
422!< jc tmp
423         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
424            ! Add asselin part
425            DO jn = 1,jpts
426               DO jk=k1,k2
427                  DO jj=j1,j2
428                     DO ji=i1,i2
429                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
430                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
431                           ztnu = tabres(ji,jj,jk,jn)
432                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
433                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
434                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
435                        ENDIF
436                     ENDDO
437                  ENDDO
438               ENDDO
439            ENDDO
440         ENDIF
441         DO jn = 1,jpts
442            DO jk=k1,k2
443               DO jj=j1,j2
444                  DO ji=i1,i2
445                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
446                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
447                     END IF
448                  END DO
449               END DO
450            END DO
451         END DO
452         !
453         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
454            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts)
455         ENDIF
456         !
457      ENDIF
458      !
459   END SUBROUTINE updateTS
460
461# endif
462
463# if defined key_vertical
464
465   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
466      !!---------------------------------------------
467      !!           *** ROUTINE updateu ***
468      !!---------------------------------------------
469      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
470      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
471      LOGICAL                                     , INTENT(in   ) :: before
472      !
473      INTEGER  ::   ji, jj, jk
474      REAL(wp) ::   zrhoy
475! VERTICAL REFINEMENT BEGIN
476      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
477      REAL(wp) :: h_in(k1:k2)
478      REAL(wp) :: h_out(1:jpk)
479      INTEGER  :: N_in, N_out
480      REAL(wp) :: h_diff, excess, thick
481      REAL(wp) :: tabin(k1:k2)
482! VERTICAL REFINEMENT END
483      !!---------------------------------------------
484      !
485      IF( before ) THEN
486         zrhoy = Agrif_Rhoy()
487         AGRIF_SpecialValue = -999._wp
488         DO jk=k1,k2
489            DO jj=j1,j2
490               DO ji=i1,i2
491                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
492                                       + (umask(ji,jj,jk)-1)*999._wp
493                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
494                                       + (umask(ji,jj,jk)-1)*999._wp
495               END DO
496            END DO
497         END DO
498      ELSE
499         tabres_child(:,:,:) = 0.
500         AGRIF_SpecialValue = 0._wp
501         DO jj=j1,j2
502            DO ji=i1,i2
503               N_in = 0
504               h_in(:) = 0._wp
505               tabin(:) = 0._wp
506               DO jk=k1,k2 !k2=jpk of child grid
507                  IF( tabres(ji,jj,jk,2) < -900) EXIT
508                  N_in = N_in + 1
509                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
510                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
511               ENDDO
512               N_out = 0
513               DO jk=1,jpk
514                  IF (umask(ji,jj,jk) == 0) EXIT
515                  N_out = N_out + 1
516                  h_out(N_out) = e3u_n(ji,jj,jk)
517               ENDDO
518               IF (N_in * N_out > 0) THEN
519                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
520                  IF (h_diff < -1.e-4) THEN
521!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
522!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
523                     excess = 0._wp
524                     DO jk=N_in,1,-1
525                        thick = MIN(-1*h_diff, h_in(jk))
526                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
527                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
528                        h_diff = h_diff + thick
529                        IF ( h_diff == 0) THEN
530                           N_in = jk
531                           h_in(jk) = h_in(jk) - thick
532                           EXIT
533                        ENDIF
534                     ENDDO
535                  ENDIF
536                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
537                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
538               ENDIF
539            ENDDO
540         ENDDO
541
542         DO jk=1,jpk
543            DO jj=j1,j2
544               DO ji=i1,i2
545                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
546                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
547                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
548                  ENDIF
549                  !
550                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
551               END DO
552            END DO
553         END DO
554      ENDIF
555      !
556   END SUBROUTINE updateu
557
558#else
559
560   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
561      !!---------------------------------------------
562      !!           *** ROUTINE updateu ***
563      !!---------------------------------------------
564      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
565      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
566      LOGICAL                                     , INTENT(in   ) :: before
567      !
568      INTEGER  :: ji, jj, jk
569      REAL(wp) :: zrhoy, zub, zunu, zuno
570      !!---------------------------------------------
571      !
572      IF( before ) THEN
573         zrhoy = Agrif_Rhoy()
574         DO jk = k1, k2
575            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
576         END DO
577      ELSE
578         DO jk=k1,k2
579            DO jj=j1,j2
580               DO ji=i1,i2
581                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 
582                  !
583                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
584                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used
585                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
586                     zunu = tabres(ji,jj,jk,1)
587                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &     
588                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
589                  ENDIF
590                  !
591                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
592               END DO
593            END DO
594         END DO
595         !
596         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
597            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
598         ENDIF
599         !
600      ENDIF
601      !
602   END SUBROUTINE updateu
603
604# endif
605
606   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
607      !!---------------------------------------------
608      !!           *** ROUTINE correct_u_bdy ***
609      !!---------------------------------------------
610      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
611      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
612      LOGICAL                                     , INTENT(in   ) :: before
613      INTEGER                                     , INTENT(in)    :: nb, ndir
614      !!
615      LOGICAL :: western_side, eastern_side 
616      !
617      INTEGER  :: jj, jk
618      REAL(wp) :: zcor
619      !!---------------------------------------------
620      !
621      IF( .NOT.before ) THEN
622         !
623         western_side  = (nb == 1).AND.(ndir == 1)
624         eastern_side  = (nb == 1).AND.(ndir == 2)
625         !
626         IF (western_side) THEN
627            DO jj=j1,j2
628               zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
629               un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
630               DO jk=1,jpkm1
631                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)
632               END DO
633            END DO
634         ENDIF
635         !
636         IF (eastern_side) THEN
637            DO jj=j1,j2
638               zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
639               un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
640               DO jk=1,jpkm1
641                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)
642               END DO
643            END DO
644         ENDIF
645         !
646      ENDIF
647      !
648   END SUBROUTINE correct_u_bdy
649
650# if  defined key_vertical
651
652   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
653      !!---------------------------------------------
654      !!           *** ROUTINE updatev ***
655      !!---------------------------------------------
656      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
657      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
658      LOGICAL                                     , INTENT(in   ) :: before
659      !
660      INTEGER  ::   ji, jj, jk
661      REAL(wp) ::   zrhox
662! VERTICAL REFINEMENT BEGIN
663      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
664      REAL(wp) :: h_in(k1:k2)
665      REAL(wp) :: h_out(1:jpk)
666      INTEGER :: N_in, N_out
667      REAL(wp) :: h_diff, excess, thick
668      REAL(wp) :: tabin(k1:k2)
669! VERTICAL REFINEMENT END
670      !!---------------------------------------------     
671      !
672      IF (before) THEN
673         zrhox = Agrif_Rhox()
674         AGRIF_SpecialValue = -999._wp
675         DO jk=k1,k2
676            DO jj=j1,j2
677               DO ji=i1,i2
678                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
679                                       + (vmask(ji,jj,jk)-1)*999._wp
680                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
681                                       + (vmask(ji,jj,jk)-1)*999._wp
682               END DO
683            END DO
684         END DO
685      ELSE
686         tabres_child(:,:,:) = 0.
687         AGRIF_SpecialValue = 0._wp
688         DO jj=j1,j2
689            DO ji=i1,i2
690               N_in = 0
691               DO jk=k1,k2
692                  IF (tabres(ji,jj,jk,2) < -900) EXIT
693                  N_in = N_in + 1
694                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
695                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
696               ENDDO
697               N_out = 0
698               DO jk=1,jpk
699                  IF (vmask(ji,jj,jk) == 0) EXIT
700                  N_out = N_out + 1
701                  h_out(N_out) = e3v_n(ji,jj,jk)
702               ENDDO
703               IF (N_in * N_out > 0) THEN
704                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
705                  IF (h_diff < -1.e-4) then
706!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
707!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
708                     excess = 0._wp
709                     DO jk=N_in,1,-1
710                        thick = MIN(-1*h_diff, h_in(jk))
711                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
712                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
713                        h_diff = h_diff + thick
714                        IF ( h_diff == 0) THEN
715                           N_in = jk
716                           h_in(jk) = h_in(jk) - thick
717                           EXIT
718                        ENDIF
719                     ENDDO
720                  ENDIF
721                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
722                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
723               ENDIF
724            ENDDO
725         ENDDO
726
727         DO jk=1,jpk
728            DO jj=j1,j2
729               DO ji=i1,i2
730                  !
731                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
732                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
733                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
734                  ENDIF
735                  !
736                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
737               END DO
738            END DO
739         END DO
740      ENDIF
741      !
742   END SUBROUTINE updatev
743
744# else
745
746   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
747      !!---------------------------------------------
748      !!           *** ROUTINE updatev ***
749      !!---------------------------------------------
750      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
751      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
752      LOGICAL                                     , INTENT(in   ) :: before
753      !
754      INTEGER  :: ji, jj, jk
755      REAL(wp) :: zrhox, zvb, zvnu, zvno
756      !!---------------------------------------------     
757      !
758      IF (before) THEN
759         zrhox = Agrif_Rhox()
760         DO jk=k1,k2
761            DO jj=j1,j2
762               DO ji=i1,i2
763                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
764               END DO
765            END DO
766         END DO
767      ELSE
768         DO jk=k1,k2
769            DO jj=j1,j2
770               DO ji=i1,i2
771                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
772                  !
773                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
774                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
775                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
776                     zvnu = tabres(ji,jj,jk,1)
777                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &     
778                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
779                  ENDIF
780                  !
781                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
782               END DO
783            END DO
784         END DO
785         !
786         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
787            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
788         ENDIF
789         !
790      ENDIF
791      !
792   END SUBROUTINE updatev
793
794# endif
795
796   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
797      !!---------------------------------------------
798      !!           *** ROUTINE correct_u_bdy ***
799      !!---------------------------------------------
800      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
801      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
802      LOGICAL                                     , INTENT(in   ) :: before
803      INTEGER                                     , INTENT(in)    :: nb, ndir
804      !!
805      LOGICAL :: southern_side, northern_side 
806      !
807      INTEGER  :: ji, jk
808      REAL(wp) :: zcor
809      !!---------------------------------------------
810      !
811      IF( .NOT.before ) THEN
812         !
813         southern_side = (nb == 2).AND.(ndir == 1)
814         northern_side = (nb == 2).AND.(ndir == 2)
815         !
816         IF (southern_side) THEN
817            DO ji=i1,i2
818               zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
819               vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
820               DO jk=1,jpkm1
821                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)
822               END DO
823            END DO
824         ENDIF
825         !
826         IF (northern_side) THEN
827            DO ji=i1,i2
828               zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
829               vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
830               DO jk=1,jpkm1
831                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)
832               END DO
833            END DO
834         ENDIF
835         !
836      ENDIF
837      !
838   END SUBROUTINE correct_v_bdy
839
840
841   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
842      !!---------------------------------------------
843      !!          *** ROUTINE updateu2d ***
844      !!---------------------------------------------
845      INTEGER, INTENT(in) :: i1, i2, j1, j2
846      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
847      LOGICAL, INTENT(in) :: before
848      !!
849      INTEGER  :: ji, jj, jk
850      REAL(wp) :: zrhoy
851      REAL(wp) :: zcorr
852      !!---------------------------------------------
853      !
854      IF (before) THEN
855         zrhoy = Agrif_Rhoy()
856         DO jj=j1,j2
857            DO ji=i1,i2
858               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
859            END DO
860         END DO
861      ELSE
862         DO jj=j1,j2
863            DO ji=i1,i2
864               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj) 
865               !   
866               ! Update "now" 3d velocities:
867               spgu(ji,jj) = 0._wp
868               DO jk=1,jpkm1
869                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
870               END DO
871               !
872               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
873               DO jk=1,jpkm1             
874                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
875               END DO
876               !
877               ! Update barotropic velocities:
878               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
879                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
880                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
881                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
882                  END IF
883               ENDIF   
884               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
885               !       
886               ! Correct "before" velocities to hold correct bt component:
887               spgu(ji,jj) = 0.e0
888               DO jk=1,jpkm1
889                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
890               END DO
891               !
892               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
893               DO jk=1,jpkm1             
894                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
895               END DO
896               !
897            END DO
898         END DO
899         !
900         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
901            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
902         ENDIF
903      ENDIF
904      !
905   END SUBROUTINE updateu2d
906
907
908   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
909      !!---------------------------------------------
910      !!          *** ROUTINE updatev2d ***
911      !!---------------------------------------------
912      INTEGER, INTENT(in) :: i1, i2, j1, j2
913      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
914      LOGICAL, INTENT(in) :: before
915      !!
916      INTEGER  :: ji, jj, jk
917      REAL(wp) :: zrhox
918      REAL(wp) :: zcorr
919      !!---------------------------------------------
920      !
921      IF (before) THEN
922         zrhox = Agrif_Rhox()
923         DO jj=j1,j2
924            DO ji=i1,i2
925               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
926            END DO
927         END DO
928      ELSE
929         DO jj=j1,j2
930            DO ji=i1,i2
931               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj) 
932               !   
933               ! Update "now" 3d velocities:
934               spgv(ji,jj) = 0.e0
935               DO jk=1,jpkm1
936                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
937               END DO
938               !
939               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
940               DO jk=1,jpkm1             
941                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
942               END DO
943               !
944               ! Update barotropic velocities:
945               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
946                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
947                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
948                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
949                  END IF
950               ENDIF             
951               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
952               !       
953               ! Correct "before" velocities to hold correct bt component:
954               spgv(ji,jj) = 0.e0
955               DO jk=1,jpkm1
956                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
957               END DO
958               !
959               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
960               DO jk=1,jpkm1             
961                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
962               END DO
963               !
964            END DO
965         END DO
966         !
967         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
968            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
969         ENDIF
970         !
971      ENDIF
972      !
973   END SUBROUTINE updatev2d
974
975
976   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
977      !!---------------------------------------------
978      !!          *** ROUTINE updateSSH ***
979      !!---------------------------------------------
980      INTEGER, INTENT(in) :: i1, i2, j1, j2
981      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
982      LOGICAL, INTENT(in) :: before
983      !!
984      INTEGER :: ji, jj
985      !!---------------------------------------------
986      !
987      IF (before) THEN
988         DO jj=j1,j2
989            DO ji=i1,i2
990               tabres(ji,jj) = sshn(ji,jj)
991            END DO
992         END DO
993      ELSE
994         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
995            DO jj=j1,j2
996               DO ji=i1,i2
997                  sshb(ji,jj) =   sshb(ji,jj) &
998                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
999               END DO
1000            END DO
1001         ENDIF
1002         !
1003         DO jj=j1,j2
1004            DO ji=i1,i2
1005               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
1006            END DO
1007         END DO
1008         !
1009         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1010            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
1011         ENDIF
1012         !
1013
1014      ENDIF
1015      !
1016   END SUBROUTINE updateSSH
1017
1018
1019   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
1020      !!---------------------------------------------
1021      !!          *** ROUTINE updateub2b ***
1022      !!---------------------------------------------
1023      INTEGER, INTENT(in) :: i1, i2, j1, j2
1024      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1025      LOGICAL, INTENT(in) :: before
1026      !!
1027      INTEGER :: ji, jj
1028      REAL(wp) :: zrhoy, za1, zcor
1029      !!---------------------------------------------
1030      !
1031      IF (before) THEN
1032         zrhoy = Agrif_Rhoy()
1033         DO jj=j1,j2
1034            DO ji=i1,i2
1035               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1036            END DO
1037         END DO
1038         tabres = zrhoy * tabres
1039      ELSE
1040         !
1041         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1042         !
1043         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1044         DO jj=j1,j2
1045            DO ji=i1,i2
1046               zcor=tabres(ji,jj) - ub2_b(ji,jj)
1047               ! Update time integrated fluxes also in case of multiply nested grids:
1048               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor 
1049               ! Update corrective fluxes:
1050               un_bf(ji,jj)  = un_bf(ji,jj) + zcor
1051               ! Update half step back fluxes:
1052               ub2_b(ji,jj) = tabres(ji,jj)
1053            END DO
1054         END DO
1055      ENDIF
1056      !
1057   END SUBROUTINE updateub2b
1058
1059   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
1060      !!---------------------------------------------
1061      !!          *** ROUTINE reflux_sshu ***
1062      !!---------------------------------------------
1063      INTEGER, INTENT(in) :: i1, i2, j1, j2
1064      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1065      LOGICAL, INTENT(in) :: before
1066      INTEGER, INTENT(in) :: nb, ndir
1067      !!
1068      LOGICAL :: western_side, eastern_side 
1069      INTEGER :: ji, jj
1070      REAL(wp) :: zrhoy, za1, zcor
1071      !!---------------------------------------------
1072      !
1073      IF (before) THEN
1074         zrhoy = Agrif_Rhoy()
1075         DO jj=j1,j2
1076            DO ji=i1,i2
1077               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1078            END DO
1079         END DO
1080         tabres = zrhoy * tabres
1081      ELSE
1082         !
1083         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1084         !
1085         western_side  = (nb == 1).AND.(ndir == 1)
1086         eastern_side  = (nb == 1).AND.(ndir == 2)
1087         !
1088         IF (western_side) THEN
1089            DO jj=j1,j2
1090               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
1091               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
1092               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor
1093            END DO
1094         ENDIF
1095         IF (eastern_side) THEN
1096            DO jj=j1,j2
1097               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
1098               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
1099               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor
1100            END DO
1101         ENDIF
1102         !
1103      ENDIF
1104      !
1105   END SUBROUTINE reflux_sshu
1106
1107   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
1108      !!---------------------------------------------
1109      !!          *** ROUTINE updatevb2b ***
1110      !!---------------------------------------------
1111      INTEGER, INTENT(in) :: i1, i2, j1, j2
1112      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1113      LOGICAL, INTENT(in) :: before
1114      !!
1115      INTEGER :: ji, jj
1116      REAL(wp) :: zrhox, za1, zcor
1117      !!---------------------------------------------
1118      !
1119      IF (before) THEN
1120         zrhox = Agrif_Rhox()
1121         DO jj=j1,j2
1122            DO ji=i1,i2
1123               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
1124            END DO
1125         END DO
1126         tabres = zrhox * tabres
1127      ELSE
1128         !
1129         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1130         !
1131         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1132         DO jj=j1,j2
1133            DO ji=i1,i2
1134               zcor=tabres(ji,jj) - vb2_b(ji,jj)
1135               ! Update time integrated fluxes also in case of multiply nested grids:
1136               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor 
1137               ! Update corrective fluxes:
1138               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor
1139               ! Update half step back fluxes:
1140               vb2_b(ji,jj) = tabres(ji,jj)
1141            END DO
1142         END DO
1143      ENDIF
1144      !
1145   END SUBROUTINE updatevb2b
1146
1147   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
1148      !!---------------------------------------------
1149      !!          *** ROUTINE reflux_sshv ***
1150      !!---------------------------------------------
1151      INTEGER, INTENT(in) :: i1, i2, j1, j2
1152      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1153      LOGICAL, INTENT(in) :: before
1154      INTEGER, INTENT(in) :: nb, ndir
1155      !!
1156      LOGICAL :: southern_side, northern_side 
1157      INTEGER :: ji, jj
1158      REAL(wp) :: zrhox, za1, zcor
1159      !!---------------------------------------------
1160      !
1161      IF (before) THEN
1162         zrhox = Agrif_Rhox()
1163         DO jj=j1,j2
1164            DO ji=i1,i2
1165               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
1166            END DO
1167         END DO
1168         tabres = zrhox * tabres
1169      ELSE
1170         !
1171         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1172         !
1173         southern_side = (nb == 2).AND.(ndir == 1)
1174         northern_side = (nb == 2).AND.(ndir == 2)
1175         !
1176         IF (southern_side) THEN
1177            DO ji=i1,i2
1178               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
1179               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
1180               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor
1181            END DO
1182         ENDIF
1183         IF (northern_side) THEN               
1184            DO ji=i1,i2
1185               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
1186               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
1187               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
1188            END DO
1189         ENDIF
1190         !
1191      ENDIF
1192      !
1193   END SUBROUTINE reflux_sshv
1194
1195   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
1196      ! currently not used
1197      !!---------------------------------------------
1198      !!           *** ROUTINE updateT ***
1199      !!---------------------------------------------
1200      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
1201      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
1202      LOGICAL, iNTENT(in) :: before
1203      !
1204      INTEGER :: ji,jj,jk
1205      REAL(wp) :: ztemp
1206      !!---------------------------------------------
1207
1208      IF (before) THEN
1209         DO jk=k1,k2
1210            DO jj=j1,j2
1211               DO ji=i1,i2
1212                  tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
1213                  tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
1214                  tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
1215               END DO
1216            END DO
1217         END DO
1218         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
1219         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
1220         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
1221      ELSE
1222         DO jk=k1,k2
1223            DO jj=j1,j2
1224               DO ji=i1,i2
1225                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
1226                     print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
1227                     print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
1228                     print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
1229                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
1230                     print *,'CORR = ',ztemp-1.
1231                     print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
1232                           tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
1233                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
1234                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
1235                  END IF
1236               END DO
1237            END DO
1238         END DO
1239      ENDIF
1240      !
1241   END SUBROUTINE update_scales
1242
1243# if defined key_zdftke
1244
1245   SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
1246      !!---------------------------------------------
1247      !!           *** ROUTINE updateen ***
1248      !!---------------------------------------------
1249      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1250      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1251      LOGICAL, INTENT(in) :: before
1252      !!---------------------------------------------
1253      !
1254      IF (before) THEN
1255         ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
1256      ELSE
1257         en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1258      ENDIF
1259      !
1260   END SUBROUTINE updateEN
1261
1262
1263   SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
1264      !!---------------------------------------------
1265      !!           *** ROUTINE updateavt ***
1266      !!---------------------------------------------
1267      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1268      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1269      LOGICAL, INTENT(in) :: before
1270      !!---------------------------------------------
1271      !
1272      IF (before) THEN
1273         ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
1274      ELSE
1275         avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1276      ENDIF
1277      !
1278   END SUBROUTINE updateAVT
1279
1280
1281   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
1282      !!---------------------------------------------
1283      !!           *** ROUTINE updateavm ***
1284      !!---------------------------------------------
1285      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1286      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1287      LOGICAL, INTENT(in) :: before
1288      !!---------------------------------------------
1289      !
1290      IF (before) THEN
1291         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1292      ELSE
1293         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1294      ENDIF
1295      !
1296   END SUBROUTINE updateAVM
1297
1298# endif /* key_zdftke */ 
1299
1300   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
1301      !!---------------------------------------------
1302      !!           *** ROUTINE updatee3t ***
1303      !!---------------------------------------------
1304      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1305      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1306      LOGICAL, INTENT(in) :: before
1307      !
1308      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1309      INTEGER :: ji,jj,jk
1310      REAL(wp) :: zcoef
1311      !!---------------------------------------------
1312      !
1313      IF (.NOT.before) THEN
1314         !
1315         ALLOCATE(ptab(i1:i2,j1:j2,1:jpk)) 
1316         !
1317         ! Update e3t from ssh (z* case only)
1318         DO jk = 1, jpkm1
1319            DO jj=j1,j2
1320               DO ji=i1,i2
1321                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
1322                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1323               END DO
1324            END DO
1325         END DO
1326         !
1327         ! 1) Updates at BEFORE time step:
1328         ! -------------------------------
1329         !
1330         ! Save "old" scale factor (prior update) for subsequent asselin correction
1331         ! of prognostic variables
1332         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
1333
1334         ! One should also save e3t_b, but lacking of workspace...
1335!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1)
1336
1337         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
1338            DO jk = 1, jpkm1
1339               DO jj=j1,j2
1340                  DO ji=i1,i2
1341                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
1342                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
1343                  END DO
1344               END DO
1345            END DO
1346            !
1347            e3w_b  (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
1348            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
1349            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
1350            !
1351            DO jk = 2, jpk
1352               DO jj = j1,j2
1353                  DO ji = i1,i2           
1354                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1355                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
1356                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
1357                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
1358                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1359                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
1360                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
1361                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
1362                  END DO
1363               END DO
1364            END DO
1365            !
1366         ENDIF       
1367         !
1368         ! 2) Updates at NOW time step:
1369         ! ----------------------------
1370         !
1371         ! Update vertical scale factor at T-points:
1372         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
1373         !
1374         ! Update total depth:
1375         ht_n(i1:i2,j1:j2) = 0._wp
1376         DO jk = 1, jpkm1
1377            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk)
1378         END DO
1379         !
1380         ! Update vertical scale factor at W-points and depths:
1381         e3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1)
1382         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
1383         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
1384         gde3w_n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh
1385         !
1386         DO jk = 2, jpk
1387            DO jj = j1,j2
1388               DO ji = i1,i2           
1389               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1390               e3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   &
1391               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1392               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
1393               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
1394                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
1395               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh
1396               END DO
1397            END DO
1398         END DO
1399         !
1400         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1401            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
1402            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
1403            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
1404            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
1405         ENDIF
1406         !
1407         DEALLOCATE(ptab)
1408      ENDIF
1409      !
1410   END SUBROUTINE updatee3t
1411
1412#else
1413CONTAINS
1414   SUBROUTINE agrif_opa_update_empty
1415      !!---------------------------------------------
1416      !!   *** ROUTINE agrif_opa_update_empty ***
1417      !!---------------------------------------------
1418      WRITE(*,*)  'agrif_opa_update : You should not have seen this print! error?'
1419   END SUBROUTINE agrif_opa_update_empty
1420#endif
1421END MODULE agrif_opa_update
1422
Note: See TracBrowser for help on using the repository browser.