source: NEMO/trunk/src/NST/agrif_oce_update.F90 @ 9748

Last change on this file since 9748 was 9598, checked in by nicolasmartin, 3 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks

  • Property svn:keywords set to Id
File size: 54.8 KB
Line 
1#define TWO_WAY        /* TWO WAY NESTING */
2#undef DECAL_FEEDBACK  /* SEPARATION of INTERFACES*/
3#undef VOL_REFLUX      /* VOLUME REFLUXING*/
4 
5MODULE agrif_oce_update
6   !!======================================================================
7   !!                   ***  MODULE  agrif_oce_interp  ***
8   !! AGRIF: update package for the ocean dynamics (OPA)
9   !!======================================================================
10   !! History :  2.0  !  2002-06  (L. Debreu)  Original code
11   !!            3.2  !  2009-04  (R. Benshila)
12   !!            3.6  !  2014-09  (R. Benshila)
13   !!----------------------------------------------------------------------
14#if defined key_agrif 
15   !!----------------------------------------------------------------------
16   !!   'key_agrif'                                              AGRIF zoom
17   !!----------------------------------------------------------------------
18   USE par_oce
19   USE oce
20   USE dom_oce
21   USE zdf_oce        ! vertical physics: ocean variables
22   USE agrif_oce
23   !
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! MPP library
26   USE domvvl         ! Need interpolation routines
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh
32   PUBLIC   Update_Scales
33
34   !!----------------------------------------------------------------------
35   !! NEMO/NST 4.0 , NEMO Consortium (2018)
36   !! $Id$
37   !! Software governed by the CeCILL licence (./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE Agrif_Update_Tra( )
42      !!----------------------------------------------------------------------
43      !!                   *** ROUTINE Agrif_Update_Tra ***
44      !!----------------------------------------------------------------------
45      !
46      IF (Agrif_Root()) RETURN
47      !
48#if defined TWO_WAY 
49      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed()
50
51      Agrif_UseSpecialValueInUpdate = .TRUE.
52      Agrif_SpecialValueFineGrid    = 0._wp
53      !
54# if ! defined DECAL_FEEDBACK
55      CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
56! near boundary update:
57!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
58# else
59      CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
60! near boundary update:
61!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
62# endif
63      !
64      Agrif_UseSpecialValueInUpdate = .FALSE.
65      !
66#endif
67      !
68   END SUBROUTINE Agrif_Update_Tra
69
70   SUBROUTINE Agrif_Update_Dyn( )
71      !!----------------------------------------------------------------------
72      !!                   *** ROUTINE Agrif_Update_Dyn ***
73      !!----------------------------------------------------------------------
74      !
75      IF (Agrif_Root()) RETURN
76      !
77#if defined TWO_WAY
78      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed()
79
80      Agrif_UseSpecialValueInUpdate = .FALSE.
81      Agrif_SpecialValueFineGrid = 0.
82      !     
83# if ! defined DECAL_FEEDBACK
84      CALL Agrif_Update_Variable(un_update_id,procname = updateU)
85      CALL Agrif_Update_Variable(vn_update_id,procname = updateV)
86! near boundary update:
87!      CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
88!      CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)
89# else
90      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
91      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
92! near boundary update:
93!      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
94!      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
95# endif
96
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
104      !
105# if ! defined DECAL_FEEDBACK
106      ! Account for updated thicknesses at boundary edges
107      IF (.NOT.ln_linssh) THEN
108!         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)
109!         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)
110      ENDIF
111# endif
112      !
113      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN
114         ! Update time integrated transports
115#  if ! defined DECAL_FEEDBACK
116         CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
117         CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
118#  else
119         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
120         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
121#  endif
122      END IF
123#endif
124      !
125   END SUBROUTINE Agrif_Update_Dyn
126
127   SUBROUTINE Agrif_Update_ssh( )
128      !!---------------------------------------------
129      !!   *** ROUTINE Agrif_Update_ssh ***
130      !!---------------------------------------------
131      !
132      IF (Agrif_Root()) RETURN
133      !
134#if defined TWO_WAY
135      !
136      Agrif_UseSpecialValueInUpdate = .TRUE.
137      Agrif_SpecialValueFineGrid = 0.
138# if ! defined DECAL_FEEDBACK
139      CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
140# else
141      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
142# endif
143      !
144      Agrif_UseSpecialValueInUpdate = .FALSE.
145      !
146#  if defined DECAL_FEEDBACK && defined VOL_REFLUX
147      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
148         ! Refluxing on ssh:
149         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)
150         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv)
151      END IF
152#  endif
153      !
154#endif
155      !
156   END SUBROUTINE Agrif_Update_ssh
157
158
159   SUBROUTINE Agrif_Update_Tke( kt )
160      !!---------------------------------------------
161      !!   *** ROUTINE Agrif_Update_Tke ***
162      !!---------------------------------------------
163      !!
164      INTEGER, INTENT(in) :: kt 
165      !
166      IF (Agrif_Root()) RETURN
167      !       
168      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
169#  if defined TWO_WAY
170
171      Agrif_UseSpecialValueInUpdate = .TRUE.
172      Agrif_SpecialValueFineGrid = 0.
173
174      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  )
175      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
176      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
177
178      Agrif_UseSpecialValueInUpdate = .FALSE.
179
180#  endif
181     
182   END SUBROUTINE Agrif_Update_Tke
183
184
185   SUBROUTINE Agrif_Update_vvl( )
186      !!---------------------------------------------
187      !!   *** ROUTINE Agrif_Update_vvl ***
188      !!---------------------------------------------
189      !
190      IF (Agrif_Root()) RETURN
191      !
192#if defined TWO_WAY 
193      !
194      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step()
195      !
196      Agrif_UseSpecialValueInUpdate = .TRUE.
197      Agrif_SpecialValueFineGrid = 0.
198      !
199      ! No interface separation here, update vertical grid at T points
200      ! everywhere over the overlapping regions (one account for refluxing in that case):
201      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t) 
202      !
203      Agrif_UseSpecialValueInUpdate = .FALSE.
204      !
205      CALL Agrif_ChildGrid_To_ParentGrid()
206      CALL dom_vvl_update_UVF
207      CALL Agrif_ParentGrid_To_ChildGrid()
208      !
209#endif
210      !
211   END SUBROUTINE Agrif_Update_vvl
212
213   SUBROUTINE dom_vvl_update_UVF
214      !!---------------------------------------------
215      !!       *** ROUTINE dom_vvl_update_UVF ***
216      !!---------------------------------------------
217      !!
218      INTEGER :: jk
219      REAL(wp):: zcoef
220      !!---------------------------------------------
221
222      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', &
223                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step()
224
225      ! Save "old" scale factor (prior update) for subsequent asselin correction
226      ! of prognostic variables
227      ! -----------------------
228      !
229      e3u_a(:,:,:) = e3u_n(:,:,:)
230      e3v_a(:,:,:) = e3v_n(:,:,:)
231!      ua(:,:,:) = e3u_b(:,:,:)
232!      va(:,:,:) = e3v_b(:,:,:)
233      hu_a(:,:) = hu_n(:,:)
234      hv_a(:,:) = hv_n(:,:)
235
236      ! 1) NOW fields
237      !--------------
238     
239         ! Vertical scale factor interpolations
240         ! ------------------------------------
241      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' )
242      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' )
243      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' )
244
245      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
246      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
247
248         ! Update total depths:
249         ! --------------------
250      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points
251      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points
252      DO jk = 1, jpkm1
253         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
254         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
255      END DO
256      !                                        ! Inverse of the local depth
257      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
258      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
259
260
261      ! 2) BEFORE fields:
262      !------------------
263      IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
264         !
265         ! Vertical scale factor interpolations
266         ! ------------------------------------
267         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  )
268         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  )
269
270         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
271         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
272
273         ! Update total depths:
274         ! --------------------
275         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points
276         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points
277         DO jk = 1, jpkm1
278            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
279            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
280         END DO
281         !                                     ! Inverse of the local depth
282         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
283         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
284      ENDIF
285      !
286   END SUBROUTINE dom_vvl_update_UVF
287
288#if defined key_vertical
289
290   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
291      !!----------------------------------------------------------------------
292      !!           *** ROUTINE updateT ***
293      !!---------------------------------------------
294      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
295      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
296      LOGICAL, INTENT(in) :: before
297      !!
298      INTEGER :: ji,jj,jk,jn
299      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
300      REAL(wp) :: h_in(k1:k2)
301      REAL(wp) :: h_out(1:jpk)
302      INTEGER  :: N_in, N_out
303      REAL(wp) :: h_diff
304      REAL(wp) :: tabin(k1:k2,n1:n2)
305      !!---------------------------------------------
306      !
307      IF (before) THEN
308         AGRIF_SpecialValue = -999._wp
309         zrho_xy = Agrif_rhox() * Agrif_rhoy() 
310         DO jn = n1,n2-1
311            DO jk=k1,k2
312               DO jj=j1,j2
313                  DO ji=i1,i2
314                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) &
315                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
316                  END DO
317               END DO
318            END DO
319         END DO
320         DO jk=k1,k2
321            DO jj=j1,j2
322               DO ji=i1,i2
323                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) &
324                                           + (tmask(ji,jj,jk)-1)*999._wp
325               END DO
326            END DO
327         END DO
328      ELSE
329         tabres_child(:,:,:,:) = 0.
330         AGRIF_SpecialValue = 0._wp
331         DO jj=j1,j2
332            DO ji=i1,i2
333               N_in = 0
334               DO jk=k1,k2 !k2 = jpk of child grid
335                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT
336                  N_in = N_in + 1
337                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
338                  h_in(N_in) = tabres(ji,jj,jk,n2)
339               ENDDO
340               N_out = 0
341               DO jk=1,jpk ! jpk of parent grid
342                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
343                  N_out = N_out + 1
344                  h_out(N_out) = e3t_n(ji,jj,jk) 
345               ENDDO
346               IF (N_in > 0) THEN !Remove this?
347                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
348                  IF (h_diff < -1.e-4) THEN
349                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
350                     print *,h_in(1:N_in)
351                     print *,h_out(1:N_out)
352                     STOP
353                  ENDIF
354                  DO jn=n1,n2-1
355                     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)
356                  ENDDO
357               ENDIF
358            ENDDO
359         ENDDO
360
361         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
362            ! Add asselin part
363            DO jn = n1,n2-1
364               DO jk=1,jpk
365                  DO jj=j1,j2
366                     DO ji=i1,i2
367                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
368                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
369                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
370                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
371                        ENDIF
372                     ENDDO
373                  ENDDO
374               ENDDO
375            ENDDO
376         ENDIF
377         DO jn = n1,n2-1
378            DO jk=1,jpk
379               DO jj=j1,j2
380                  DO ji=i1,i2
381                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
382                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
383                     END IF
384                  END DO
385               END DO
386            END DO
387         END DO
388      ENDIF
389      !
390   END SUBROUTINE updateTS
391
392# else
393
394   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
395      !!---------------------------------------------
396      !!           *** ROUTINE updateT ***
397      !!---------------------------------------------
398      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
399      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
400      LOGICAL, INTENT(in) :: before
401      !!
402      INTEGER :: ji,jj,jk,jn
403      REAL(wp) :: ztb, ztnu, ztno
404      !!---------------------------------------------
405      !
406      IF (before) THEN
407         DO jn = 1,jpts
408            DO jk=k1,k2
409               DO jj=j1,j2
410                  DO ji=i1,i2
411!> jc tmp
412                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk)
413!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk)
414!< jc tmp
415                  END DO
416               END DO
417            END DO
418         END DO
419      ELSE
420!> jc tmp
421         DO jn = 1,jpts
422            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
423                                         & * tmask(i1:i2,j1:j2,k1:k2)
424         ENDDO
425!< jc tmp
426         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
427            ! Add asselin part
428            DO jn = 1,jpts
429               DO jk = k1, k2
430                  DO jj = j1, j2
431                     DO ji = i1, i2
432                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
433                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used
434                           ztnu = tabres(ji,jj,jk,jn)
435                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk)
436                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  & 
437                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk)
438                        ENDIF
439                     END DO
440                  END DO
441               END DO
442            END DO
443         ENDIF
444         DO jn = 1,jpts
445            DO jk=k1,k2
446               DO jj=j1,j2
447                  DO ji=i1,i2
448                     IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
449                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk)
450                     END IF
451                  END DO
452               END DO
453            END DO
454         END DO
455         !
456         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
457            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts)
458         ENDIF
459         !
460      ENDIF
461      !
462   END SUBROUTINE updateTS
463
464# endif
465
466# if defined key_vertical
467
468   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
469      !!---------------------------------------------
470      !!           *** ROUTINE updateu ***
471      !!---------------------------------------------
472      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
473      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
474      LOGICAL                                     , INTENT(in   ) :: before
475      !
476      INTEGER ::   ji, jj, jk
477      REAL(wp)::   zrhoy
478! VERTICAL REFINEMENT BEGIN
479      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
480      REAL(wp) :: h_in(k1:k2)
481      REAL(wp) :: h_out(1:jpk)
482      INTEGER  :: N_in, N_out
483      REAL(wp) :: h_diff, excess, thick
484      REAL(wp) :: tabin(k1:k2)
485! VERTICAL REFINEMENT END
486      !!---------------------------------------------
487      !
488      IF( before ) THEN
489         zrhoy = Agrif_Rhoy()
490         AGRIF_SpecialValue = -999._wp
491         DO jk=k1,k2
492            DO jj=j1,j2
493               DO ji=i1,i2
494                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
495                                       + (umask(ji,jj,jk)-1)*999._wp
496                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
497                                       + (umask(ji,jj,jk)-1)*999._wp
498               END DO
499            END DO
500         END DO
501      ELSE
502         tabres_child(:,:,:) = 0.
503         AGRIF_SpecialValue = 0._wp
504         DO jj=j1,j2
505            DO ji=i1,i2
506               N_in = 0
507               h_in(:) = 0._wp
508               tabin(:) = 0._wp
509               DO jk=k1,k2 !k2=jpk of child grid
510                  IF( tabres(ji,jj,jk,2) < -900) EXIT
511                  N_in = N_in + 1
512                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
513                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
514               ENDDO
515               N_out = 0
516               DO jk=1,jpk
517                  IF (umask(ji,jj,jk) == 0) EXIT
518                  N_out = N_out + 1
519                  h_out(N_out) = e3u_n(ji,jj,jk)
520               ENDDO
521               IF (N_in * N_out > 0) THEN
522                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
523                  IF (h_diff < -1.e-4) THEN
524!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
525!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
526                     excess = 0._wp
527                     DO jk=N_in,1,-1
528                        thick = MIN(-1*h_diff, h_in(jk))
529                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
530                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
531                        h_diff = h_diff + thick
532                        IF ( h_diff == 0) THEN
533                           N_in = jk
534                           h_in(jk) = h_in(jk) - thick
535                           EXIT
536                        ENDIF
537                     ENDDO
538                  ENDIF
539                  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)
540                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
541               ENDIF
542            ENDDO
543         ENDDO
544
545         DO jk=1,jpk
546            DO jj=j1,j2
547               DO ji=i1,i2
548                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
549                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
550                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
551                  ENDIF
552                  !
553                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
554               END DO
555            END DO
556         END DO
557      ENDIF
558      !
559   END SUBROUTINE updateu
560
561#else
562
563   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
564      !!---------------------------------------------
565      !!           *** ROUTINE updateu ***
566      !!---------------------------------------------
567      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
568      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
569      LOGICAL                                     , INTENT(in   ) :: before
570      !
571      INTEGER  :: ji, jj, jk
572      REAL(wp) :: zrhoy, zub, zunu, zuno
573      !!---------------------------------------------
574      !
575      IF( before ) THEN
576         zrhoy = Agrif_Rhoy()
577         DO jk = k1, k2
578            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)
579         END DO
580      ELSE
581         DO jk=k1,k2
582            DO jj=j1,j2
583               DO ji=i1,i2
584                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 
585                  !
586                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
587                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used
588                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk)
589                     zunu = tabres(ji,jj,jk,1)
590                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &     
591                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk)
592                  ENDIF
593                  !
594                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)
595               END DO
596            END DO
597         END DO
598         !
599         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
600            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2)
601         ENDIF
602         !
603      ENDIF
604      !
605   END SUBROUTINE updateu
606
607# endif
608
609   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
610      !!---------------------------------------------
611      !!           *** ROUTINE correct_u_bdy ***
612      !!---------------------------------------------
613      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
614      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
615      LOGICAL                                     , INTENT(in   ) :: before
616      INTEGER                                     , INTENT(in)    :: nb, ndir
617      !!
618      LOGICAL :: western_side, eastern_side 
619      !
620      INTEGER  :: jj, jk
621      REAL(wp) :: zcor
622      !!---------------------------------------------
623      !
624      IF( .NOT.before ) THEN
625         !
626         western_side  = (nb == 1).AND.(ndir == 1)
627         eastern_side  = (nb == 1).AND.(ndir == 2)
628         !
629         IF (western_side) THEN
630            DO jj=j1,j2
631               zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj)
632               un_b(i1-1,jj) = un_b(i1-1,jj) + zcor
633               DO jk=1,jpkm1
634                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk)
635               END DO
636            END DO
637         ENDIF
638         !
639         IF (eastern_side) THEN
640            DO jj=j1,j2
641               zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj)
642               un_b(i2+1,jj) = un_b(i2+1,jj) + zcor
643               DO jk=1,jpkm1
644                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk)
645               END DO
646            END DO
647         ENDIF
648         !
649      ENDIF
650      !
651   END SUBROUTINE correct_u_bdy
652
653# if  defined key_vertical
654
655   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
656      !!---------------------------------------------
657      !!           *** ROUTINE updatev ***
658      !!---------------------------------------------
659      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
660      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
661      LOGICAL                                     , INTENT(in   ) :: before
662      !
663      INTEGER  ::   ji, jj, jk
664      REAL(wp) ::   zrhox
665! VERTICAL REFINEMENT BEGIN
666      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
667      REAL(wp) :: h_in(k1:k2)
668      REAL(wp) :: h_out(1:jpk)
669      INTEGER :: N_in, N_out
670      REAL(wp) :: h_diff, excess, thick
671      REAL(wp) :: tabin(k1:k2)
672! VERTICAL REFINEMENT END
673      !!---------------------------------------------     
674      !
675      IF( before ) THEN
676         zrhox = Agrif_Rhox()
677         AGRIF_SpecialValue = -999._wp
678         DO jk=k1,k2
679            DO jj=j1,j2
680               DO ji=i1,i2
681                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
682                                       + (vmask(ji,jj,jk)-1)*999._wp
683                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
684                                       + (vmask(ji,jj,jk)-1)*999._wp
685               END DO
686            END DO
687         END DO
688      ELSE
689         tabres_child(:,:,:) = 0.
690         AGRIF_SpecialValue = 0._wp
691         DO jj=j1,j2
692            DO ji=i1,i2
693               N_in = 0
694               DO jk=k1,k2
695                  IF (tabres(ji,jj,jk,2) < -900) EXIT
696                  N_in = N_in + 1
697                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
698                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
699               ENDDO
700               N_out = 0
701               DO jk=1,jpk
702                  IF (vmask(ji,jj,jk) == 0) EXIT
703                  N_out = N_out + 1
704                  h_out(N_out) = e3v_n(ji,jj,jk)
705               ENDDO
706               IF (N_in * N_out > 0) THEN
707                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
708                  IF (h_diff < -1.e-4) then
709!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
710!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
711                     excess = 0._wp
712                     DO jk=N_in,1,-1
713                        thick = MIN(-1*h_diff, h_in(jk))
714                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
715                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
716                        h_diff = h_diff + thick
717                        IF ( h_diff == 0) THEN
718                           N_in = jk
719                           h_in(jk) = h_in(jk) - thick
720                           EXIT
721                        ENDIF
722                     ENDDO
723                  ENDIF
724                  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)
725                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
726               ENDIF
727            ENDDO
728         ENDDO
729
730         DO jk=1,jpk
731            DO jj=j1,j2
732               DO ji=i1,i2
733                  !
734                  IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part
735                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
736                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
737                  ENDIF
738                  !
739                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
740               END DO
741            END DO
742         END DO
743      ENDIF
744      !
745   END SUBROUTINE updatev
746
747# else
748
749   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
750      !!---------------------------------------------
751      !!           *** ROUTINE updatev ***
752      !!---------------------------------------------
753      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
754      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
755      LOGICAL                                     , INTENT(in   ) :: before
756      !
757      INTEGER  :: ji, jj, jk
758      REAL(wp) :: zrhox, zvb, zvnu, zvno
759      !!---------------------------------------------     
760      !
761      IF (before) THEN
762         zrhox = Agrif_Rhox()
763         DO jk=k1,k2
764            DO jj=j1,j2
765               DO ji=i1,i2
766                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
767               END DO
768            END DO
769         END DO
770      ELSE
771         DO jk=k1,k2
772            DO jj=j1,j2
773               DO ji=i1,i2
774                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
775                  !
776                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
777                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used
778                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk)
779                     zvnu = tabres(ji,jj,jk,1)
780                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &     
781                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk)
782                  ENDIF
783                  !
784                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)
785               END DO
786            END DO
787         END DO
788         !
789         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
790            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2)
791         ENDIF
792         !
793      ENDIF
794      !
795   END SUBROUTINE updatev
796
797# endif
798
799   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
800      !!---------------------------------------------
801      !!           *** ROUTINE correct_u_bdy ***
802      !!---------------------------------------------
803      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
804      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
805      LOGICAL                                     , INTENT(in   ) :: before
806      INTEGER                                     , INTENT(in)    :: nb, ndir
807      !!
808      LOGICAL :: southern_side, northern_side 
809      !
810      INTEGER  :: ji, jk
811      REAL(wp) :: zcor
812      !!---------------------------------------------
813      !
814      IF( .NOT.before ) THEN
815         !
816         southern_side = (nb == 2).AND.(ndir == 1)
817         northern_side = (nb == 2).AND.(ndir == 2)
818         !
819         IF (southern_side) THEN
820            DO ji=i1,i2
821               zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1)
822               vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor
823               DO jk=1,jpkm1
824                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk)
825               END DO
826            END DO
827         ENDIF
828         !
829         IF (northern_side) THEN
830            DO ji=i1,i2
831               zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1)
832               vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor
833               DO jk=1,jpkm1
834                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk)
835               END DO
836            END DO
837         ENDIF
838         !
839      ENDIF
840      !
841   END SUBROUTINE correct_v_bdy
842
843
844   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
845      !!----------------------------------------------------------------------
846      !!                      *** ROUTINE updateu2d ***
847      !!----------------------------------------------------------------------
848      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
849      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
850      LOGICAL                         , INTENT(in   ) ::   before
851      !!
852      INTEGER  :: ji, jj, jk
853      REAL(wp) :: zrhoy
854      REAL(wp) :: zcorr
855      !!---------------------------------------------
856      !
857      IF( before ) THEN
858         zrhoy = Agrif_Rhoy()
859         DO jj=j1,j2
860            DO ji=i1,i2
861               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
862            END DO
863         END DO
864      ELSE
865         DO jj=j1,j2
866            DO ji=i1,i2
867               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj) 
868               !   
869               ! Update "now" 3d velocities:
870               spgu(ji,jj) = 0._wp
871               DO jk=1,jpkm1
872                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
873               END DO
874               !
875               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj)
876               DO jk=1,jpkm1             
877                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
878               END DO
879               !
880               ! Update barotropic velocities:
881               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
882                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
883                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj)
884                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
885                  END IF
886               ENDIF   
887               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1)
888               !       
889               ! Correct "before" velocities to hold correct bt component:
890               spgu(ji,jj) = 0.e0
891               DO jk=1,jpkm1
892                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
893               END DO
894               !
895               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj)
896               DO jk=1,jpkm1             
897                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
898               END DO
899               !
900            END DO
901         END DO
902         !
903         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
904            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2)
905         ENDIF
906      ENDIF
907      !
908   END SUBROUTINE updateu2d
909
910
911   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
912      !!----------------------------------------------------------------------
913      !!                   *** ROUTINE updatev2d ***
914      !!----------------------------------------------------------------------
915      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
916      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
917      LOGICAL                         , INTENT(in   ) ::   before
918      !
919      INTEGER  :: ji, jj, jk
920      REAL(wp) :: zrhox, zcorr
921      !!----------------------------------------------------------------------
922      !
923      IF( before ) THEN
924         zrhox = Agrif_Rhox()
925         DO jj=j1,j2
926            DO ji=i1,i2
927               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
928            END DO
929         END DO
930      ELSE
931         DO jj=j1,j2
932            DO ji=i1,i2
933               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj) 
934               !   
935               ! Update "now" 3d velocities:
936               spgv(ji,jj) = 0.e0
937               DO jk=1,jpkm1
938                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
939               END DO
940               !
941               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj)
942               DO jk=1,jpkm1             
943                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
944               END DO
945               !
946               ! Update barotropic velocities:
947               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
948                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
949                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj)
950                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
951                  END IF
952               ENDIF             
953               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1)
954               !       
955               ! Correct "before" velocities to hold correct bt component:
956               spgv(ji,jj) = 0.e0
957               DO jk=1,jpkm1
958                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
959               END DO
960               !
961               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj)
962               DO jk=1,jpkm1             
963                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
964               END DO
965               !
966            END DO
967         END DO
968         !
969         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
970            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2)
971         ENDIF
972         !
973      ENDIF
974      !
975   END SUBROUTINE updatev2d
976
977
978   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
979      !!----------------------------------------------------------------------
980      !!                   *** ROUTINE updateSSH ***
981      !!----------------------------------------------------------------------
982      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
983      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
984      LOGICAL                         , INTENT(in   ) ::   before
985      !!
986      INTEGER :: ji, jj
987      !!----------------------------------------------------------------------
988      !
989      IF( before ) THEN
990         DO jj=j1,j2
991            DO ji=i1,i2
992               tabres(ji,jj) = sshn(ji,jj)
993            END DO
994         END DO
995      ELSE
996         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
997            DO jj=j1,j2
998               DO ji=i1,i2
999                  sshb(ji,jj) =   sshb(ji,jj) &
1000                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
1001               END DO
1002            END DO
1003         ENDIF
1004         !
1005         DO jj=j1,j2
1006            DO ji=i1,i2
1007               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
1008            END DO
1009         END DO
1010         !
1011         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1012            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2)
1013         ENDIF
1014         !
1015
1016      ENDIF
1017      !
1018   END SUBROUTINE updateSSH
1019
1020
1021   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
1022      !!----------------------------------------------------------------------
1023      !!                      *** ROUTINE updateub2b ***
1024      !!----------------------------------------------------------------------
1025      INTEGER                            , INTENT(in) ::   i1, i2, j1, j2
1026      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
1027      LOGICAL                            , INTENT(in) ::   before
1028      !!
1029      INTEGER :: ji, jj
1030      REAL(wp) :: zrhoy, za1, zcor
1031      !!---------------------------------------------
1032      !
1033      IF (before) THEN
1034         zrhoy = Agrif_Rhoy()
1035         DO jj=j1,j2
1036            DO ji=i1,i2
1037               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1038            END DO
1039         END DO
1040         tabres = zrhoy * tabres
1041      ELSE
1042         !
1043         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1044         !
1045         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1046         DO jj=j1,j2
1047            DO ji=i1,i2
1048               zcor=tabres(ji,jj) - ub2_b(ji,jj)
1049               ! Update time integrated fluxes also in case of multiply nested grids:
1050               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor 
1051               ! Update corrective fluxes:
1052               un_bf(ji,jj)  = un_bf(ji,jj) + zcor
1053               ! Update half step back fluxes:
1054               ub2_b(ji,jj) = tabres(ji,jj)
1055            END DO
1056         END DO
1057      ENDIF
1058      !
1059   END SUBROUTINE updateub2b
1060
1061   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
1062      !!---------------------------------------------
1063      !!          *** ROUTINE reflux_sshu ***
1064      !!---------------------------------------------
1065      INTEGER, INTENT(in) :: i1, i2, j1, j2
1066      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1067      LOGICAL, INTENT(in) :: before
1068      INTEGER, INTENT(in) :: nb, ndir
1069      !!
1070      LOGICAL :: western_side, eastern_side 
1071      INTEGER :: ji, jj
1072      REAL(wp) :: zrhoy, za1, zcor
1073      !!---------------------------------------------
1074      !
1075      IF (before) THEN
1076         zrhoy = Agrif_Rhoy()
1077         DO jj=j1,j2
1078            DO ji=i1,i2
1079               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
1080            END DO
1081         END DO
1082         tabres = zrhoy * tabres
1083      ELSE
1084         !
1085         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
1086         !
1087         western_side  = (nb == 1).AND.(ndir == 1)
1088         eastern_side  = (nb == 1).AND.(ndir == 2)
1089         !
1090         IF (western_side) THEN
1091            DO jj=j1,j2
1092               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
1093               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
1094               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor
1095            END DO
1096         ENDIF
1097         IF (eastern_side) THEN
1098            DO jj=j1,j2
1099               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
1100               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
1101               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor
1102            END DO
1103         ENDIF
1104         !
1105      ENDIF
1106      !
1107   END SUBROUTINE reflux_sshu
1108
1109   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
1110      !!----------------------------------------------------------------------
1111      !!                      *** ROUTINE updatevb2b ***
1112      !!----------------------------------------------------------------------
1113      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
1114      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
1115      LOGICAL                         , INTENT(in   ) ::   before
1116      !!
1117      INTEGER :: ji, jj
1118      REAL(wp) :: zrhox, za1, zcor
1119      !!---------------------------------------------
1120      !
1121      IF( before ) THEN
1122         zrhox = Agrif_Rhox()
1123         DO jj=j1,j2
1124            DO ji=i1,i2
1125               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
1126            END DO
1127         END DO
1128         tabres = zrhox * tabres
1129      ELSE
1130         !
1131         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1132         !
1133         za1 = 1._wp / REAL(Agrif_rhot(), wp)
1134         DO jj=j1,j2
1135            DO ji=i1,i2
1136               zcor=tabres(ji,jj) - vb2_b(ji,jj)
1137               ! Update time integrated fluxes also in case of multiply nested grids:
1138               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor 
1139               ! Update corrective fluxes:
1140               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor
1141               ! Update half step back fluxes:
1142               vb2_b(ji,jj) = tabres(ji,jj)
1143            END DO
1144         END DO
1145      ENDIF
1146      !
1147   END SUBROUTINE updatevb2b
1148
1149   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
1150      !!---------------------------------------------
1151      !!          *** ROUTINE reflux_sshv ***
1152      !!---------------------------------------------
1153      INTEGER, INTENT(in) :: i1, i2, j1, j2
1154      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
1155      LOGICAL, INTENT(in) :: before
1156      INTEGER, INTENT(in) :: nb, ndir
1157      !!
1158      LOGICAL :: southern_side, northern_side 
1159      INTEGER :: ji, jj
1160      REAL(wp) :: zrhox, za1, zcor
1161      !!---------------------------------------------
1162      !
1163      IF (before) THEN
1164         zrhox = Agrif_Rhox()
1165         DO jj=j1,j2
1166            DO ji=i1,i2
1167               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
1168            END DO
1169         END DO
1170         tabres = zrhox * tabres
1171      ELSE
1172         !
1173         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
1174         !
1175         southern_side = (nb == 2).AND.(ndir == 1)
1176         northern_side = (nb == 2).AND.(ndir == 2)
1177         !
1178         IF (southern_side) THEN
1179            DO ji=i1,i2
1180               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
1181               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
1182               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor
1183            END DO
1184         ENDIF
1185         IF (northern_side) THEN               
1186            DO ji=i1,i2
1187               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
1188               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
1189               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
1190            END DO
1191         ENDIF
1192         !
1193      ENDIF
1194      !
1195   END SUBROUTINE reflux_sshv
1196
1197   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
1198      !
1199      ! ====>>>>>>>>>>    currently not used
1200      !
1201      !!----------------------------------------------------------------------
1202      !!                      *** ROUTINE updateT ***
1203      !!----------------------------------------------------------------------
1204      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
1205      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
1206      LOGICAL                                    , INTENT(in   ) ::   before
1207      !!
1208      INTEGER :: ji,jj,jk
1209      REAL(wp) :: ztemp
1210      !!----------------------------------------------------------------------
1211
1212      IF (before) THEN
1213         DO jk=k1,k2
1214            DO jj=j1,j2
1215               DO ji=i1,i2
1216                  tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
1217                  tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
1218                  tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
1219               END DO
1220            END DO
1221         END DO
1222         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
1223         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
1224         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
1225      ELSE
1226         DO jk=k1,k2
1227            DO jj=j1,j2
1228               DO ji=i1,i2
1229                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
1230                     print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
1231                     print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
1232                     print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
1233                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
1234                     print *,'CORR = ',ztemp-1.
1235                     print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
1236                           tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
1237                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
1238                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
1239                  END IF
1240               END DO
1241            END DO
1242         END DO
1243      ENDIF
1244      !
1245   END SUBROUTINE update_scales
1246
1247
1248   SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
1249      !!----------------------------------------------------------------------
1250      !!                      *** ROUTINE updateen ***
1251      !!----------------------------------------------------------------------
1252      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1253      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1254      LOGICAL                               , INTENT(in   ) ::   before
1255      !!----------------------------------------------------------------------
1256      !
1257      IF( before ) THEN
1258         ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
1259      ELSE
1260         en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1261      ENDIF
1262      !
1263   END SUBROUTINE updateEN
1264
1265
1266   SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
1267      !!----------------------------------------------------------------------
1268      !!                      *** ROUTINE updateavt ***
1269      !!----------------------------------------------------------------------
1270      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1271      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1272      LOGICAL                               , INTENT(in   ) ::   before
1273      !!----------------------------------------------------------------------
1274      !
1275      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
1276      ELSE                ;   avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1277      ENDIF
1278      !
1279   END SUBROUTINE updateAVT
1280
1281
1282   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
1283      !!---------------------------------------------
1284      !!           *** ROUTINE updateavm ***
1285      !!----------------------------------------------------------------------
1286      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
1287      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
1288      LOGICAL                               , INTENT(in   ) ::   before
1289      !!----------------------------------------------------------------------
1290      !
1291      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1292      ELSE                ;   avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
1293      ENDIF
1294      !
1295   END SUBROUTINE updateAVM
1296
1297   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
1298      !!---------------------------------------------
1299      !!           *** ROUTINE updatee3t ***
1300      !!---------------------------------------------
1301      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
1302      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
1303      LOGICAL, INTENT(in) :: before
1304      !
1305      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
1306      INTEGER :: ji,jj,jk
1307      REAL(wp) :: zcoef
1308      !!---------------------------------------------
1309      !
1310      IF (.NOT.before) THEN
1311         !
1312         ALLOCATE(ptab(i1:i2,j1:j2,1:jpk)) 
1313         !
1314         ! Update e3t from ssh (z* case only)
1315         DO jk = 1, jpkm1
1316            DO jj=j1,j2
1317               DO ji=i1,i2
1318                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
1319                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
1320               END DO
1321            END DO
1322         END DO
1323         !
1324         ! 1) Updates at BEFORE time step:
1325         ! -------------------------------
1326         !
1327         ! Save "old" scale factor (prior update) for subsequent asselin correction
1328         ! of prognostic variables
1329         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
1330
1331         ! One should also save e3t_b, but lacking of workspace...
1332!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1)
1333
1334         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
1335            DO jk = 1, jpkm1
1336               DO jj=j1,j2
1337                  DO ji=i1,i2
1338                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) &
1339                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) )
1340                  END DO
1341               END DO
1342            END DO
1343            !
1344            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)
1345            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp
1346            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1)
1347            !
1348            DO jk = 2, jpk
1349               DO jj = j1,j2
1350                  DO ji = i1,i2           
1351                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1352                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        & 
1353                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  &
1354                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        &
1355                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1356                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
1357                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
1358                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
1359                  END DO
1360               END DO
1361            END DO
1362            !
1363         ENDIF       
1364         !
1365         ! 2) Updates at NOW time step:
1366         ! ----------------------------
1367         !
1368         ! Update vertical scale factor at T-points:
1369         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
1370         !
1371         ! Update total depth:
1372         ht_n(i1:i2,j1:j2) = 0._wp
1373         DO jk = 1, jpkm1
1374            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)
1375         END DO
1376         !
1377         ! Update vertical scale factor at W-points and depths:
1378         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)
1379         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1)
1380         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp
1381         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
1382         !
1383         DO jk = 2, jpk
1384            DO jj = j1,j2
1385               DO ji = i1,i2           
1386               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
1387               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) )   &
1388               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) )
1389               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
1390               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
1391                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
1392               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
1393               END DO
1394            END DO
1395         END DO
1396         !
1397         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
1398            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk)
1399            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk)
1400            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk)
1401            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk)
1402         ENDIF
1403         !
1404         DEALLOCATE(ptab)
1405      ENDIF
1406      !
1407   END SUBROUTINE updatee3t
1408
1409#else
1410   !!----------------------------------------------------------------------
1411   !!   Empty module                                          no AGRIF zoom
1412   !!----------------------------------------------------------------------
1413CONTAINS
1414   SUBROUTINE agrif_oce_update_empty
1415      WRITE(*,*)  'agrif_oce_update : You should not have seen this print! error?'
1416   END SUBROUTINE agrif_oce_update_empty
1417#endif
1418
1419   !!======================================================================
1420END MODULE agrif_oce_update
1421
Note: See TracBrowser for help on using the repository browser.