New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_opa_update.F90 in branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 8690

Last change on this file since 8690 was 8690, checked in by timgraham, 7 years ago

Bug fixes to get full bathy overflow test working

  • Property svn:keywords set to Id
File size: 30.9 KB
Line 
1#define TWO_WAY        /* TWO WAY NESTING */
2#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/
3 
4MODULE agrif_opa_update
5#if defined key_agrif 
6   USE par_oce
7   USE oce
8   USE dom_oce
9   USE agrif_oce
10   USE in_out_manager  ! I/O manager
11   USE lib_mpp
12   USE wrk_nemo 
13   USE zdf_oce        ! vertical physics: ocean variables
14
15   IMPLICIT NONE
16   PRIVATE
17
18   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn,Update_Scales
19# if defined key_zdftke
20   PUBLIC Agrif_Update_Tke
21# endif
22   !!----------------------------------------------------------------------
23   !! NEMO/NST 3.6 , NEMO Consortium (2010)
24   !! $Id$
25   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
26   !!----------------------------------------------------------------------
27CONTAINS
28
29   RECURSIVE SUBROUTINE Agrif_Update_Tra( )
30      !!---------------------------------------------
31      !!   *** ROUTINE Agrif_Update_Tra ***
32      !!---------------------------------------------
33      !
34      IF (Agrif_Root()) RETURN
35      !
36#if defined TWO_WAY 
37      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed(), 'nbcline', nbcline
38
39      Agrif_UseSpecialValueInUpdate = .TRUE.
40      Agrif_SpecialValueFineGrid = 0.
41      !
42      IF (MOD(nbcline,nbclineupdate) == 0) THEN
43# if ! defined DECAL_FEEDBACK
44         CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
45# else
46         CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
47# endif
48      ELSE
49# if ! defined DECAL_FEEDBACK
50         CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
51# else
52         CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
53# endif
54      ENDIF
55      !
56      Agrif_UseSpecialValueInUpdate = .FALSE.
57      !
58      IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:
59         CALL Agrif_ChildGrid_To_ParentGrid()
60         CALL Agrif_Update_Tra()
61         CALL Agrif_ParentGrid_To_ChildGrid()
62      ENDIF
63      !
64#endif
65      !
66   END SUBROUTINE Agrif_Update_Tra
67
68
69   RECURSIVE SUBROUTINE Agrif_Update_Dyn( )
70      !!---------------------------------------------
71      !!   *** ROUTINE Agrif_Update_Dyn ***
72      !!---------------------------------------------
73      !
74      IF (Agrif_Root()) RETURN
75      !
76#if defined TWO_WAY
77      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed(), 'nbcline', nbcline
78
79      Agrif_UseSpecialValueInUpdate = .FALSE.
80      Agrif_SpecialValueFineGrid = 0.
81      !     
82      IF (mod(nbcline,nbclineupdate) == 0) THEN
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# else
87         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU)
88         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV)
89# endif
90      ELSE
91# if ! defined DECAL_FEEDBACK
92         CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU)
93         CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)         
94# else
95         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU)
96         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
97# endif
98      ENDIF
99
100# if ! defined DECAL_FEEDBACK
101      CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
102      CALL Agrif_Update_Variable(e2v_id,procname = updateV2d) 
103# else
104      CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
105      CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d) 
106# endif
107
108      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
109         ! Update time integrated transports
110         IF (mod(nbcline,nbclineupdate) == 0) THEN
111#  if ! defined DECAL_FEEDBACK
112            CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
113            CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
114#  else
115            CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
116            CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
117#  endif
118         ELSE
119#  if ! defined DECAL_FEEDBACK
120            CALL Agrif_Update_Variable(ub2b_update_id,locupdate=(/0,1/),procname = updateub2b)
121            CALL Agrif_Update_Variable(vb2b_update_id,locupdate=(/0,1/),procname = updatevb2b)
122#  else
123            CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateub2b)
124            CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updatevb2b)
125#  endif
126         ENDIF
127      END IF
128      !
129      nbcline = nbcline + 1
130      !
131      Agrif_UseSpecialValueInUpdate = .TRUE.
132      Agrif_SpecialValueFineGrid = 0.
133# if ! defined DECAL_FEEDBACK
134      CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
135# else
136      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
137# endif
138      Agrif_UseSpecialValueInUpdate = .FALSE.
139      !
140#endif
141      !
142      ! Do recursive update:
143      IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update:
144         CALL Agrif_ChildGrid_To_ParentGrid()
145         CALL Agrif_Update_Dyn()
146         CALL Agrif_ParentGrid_To_ChildGrid()
147      ENDIF
148      !
149   END SUBROUTINE Agrif_Update_Dyn
150
151# if defined key_zdftke
152
153   SUBROUTINE Agrif_Update_Tke( kt )
154      !!---------------------------------------------
155      !!   *** ROUTINE Agrif_Update_Tke ***
156      !!---------------------------------------------
157      !!
158      INTEGER, INTENT(in) :: kt
159      !       
160      return
161     
162      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
163#  if defined TWO_WAY
164
165      Agrif_UseSpecialValueInUpdate = .TRUE.
166      Agrif_SpecialValueFineGrid = 0.
167
168      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  )
169      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
170      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
171
172      Agrif_UseSpecialValueInUpdate = .FALSE.
173
174#  endif
175     
176   END SUBROUTINE Agrif_Update_Tke
177   
178# endif /* key_zdftke */
179
180   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
181      !!---------------------------------------------
182      !!           *** ROUTINE updateT ***
183      !!---------------------------------------------
184      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
185      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
186      LOGICAL, INTENT(in) :: before
187      !!
188      INTEGER :: ji,jj,jk,jn
189! VERTICAL REFINEMENT BEGIN
190      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
191      REAL(wp) :: h_in(k1:k2)
192      REAL(wp) :: h_out(1:jpk)
193      INTEGER  :: N_in, N_out
194      REAL(wp) :: h_diff
195      REAL(wp) :: zrho_xy
196      REAL(wp) :: tabin(k1:k2,n1:n2)
197! VERTICAL REFINEMENT END
198      !!---------------------------------------------
199      !
200      IF (before) THEN
201# if defined key_vertical
202            AGRIF_SpecialValue = -999._wp
203            zrho_xy = Agrif_rhox() * Agrif_rhoy() 
204            DO jn = n1,n2-1
205               DO jk=k1,k2
206                  DO jj=j1,j2
207                     DO ji=i1,i2
208                        tabres(ji,jj,jk,jn) = (zrho_xy * tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) &
209                                              * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
210                     END DO
211                  END DO
212               END DO
213            END DO
214            DO jk=k1,k2
215               DO jj=j1,j2
216                  DO ji=i1,i2
217                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)  &
218                                              + (tmask(ji,jj,jk)-1)*999._wp
219                  END DO
220               END DO
221            END DO
222#else
223            DO jn = n1,n2-1
224               DO jk=k1,k2
225                  DO jj=j1,j2
226                     DO ji=i1,i2
227                        tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
228                     END DO
229                  END DO
230               END DO
231            END DO       
232#endif
233      ELSE
234! VERTICAL REFINEMENT BEGIN
235         tabres_child(:,:,:,:) = 0.
236# if defined key_vertical
237         AGRIF_SpecialValue = 0._wp
238         DO jj=j1,j2
239         DO ji=i1,i2
240           N_in = 0
241           DO jk=k1,k2 !k2 = jpk of child grid
242             IF (tabres(ji,jj,jk,n2) > -900) EXIT
243             N_in = N_in + 1
244             tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
245             h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj)
246           ENDDO
247           N_out = 0
248           DO jk=1,jpk ! jpk of parent grid
249             IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF
250             N_out = N_out + 1
251             h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
252           ENDDO
253!           IF(ji.EQ.i1 .AND. jj.EQ.j1) print *,'1st parent point',sum(h_in(1:N_in)), sum(h_out(1:N_out))
254           IF (N_in > 0) THEN
255             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
256! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
257!             IF abs(h_diff > 1.e-8) THEN
258!               N_in = N_in + 1
259!               h_in(N_in) = h_diff
260!               tabin(N_in,:) = tabin(N_in-1,:)
261             IF (h_diff < -1.e-4) THEN
262             print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,e1e2t(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out
263             print *,h_in(1:N_in)
264             print *,h_out(1:N_out)
265             STOP
266!               N_out = N_out + 1
267!               h_out(N_out) = - h_diff
268             ENDIF
269             DO jn=n1,n2-1
270               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)
271             ENDDO
272          ENDIF
273         ENDDO
274         ENDDO
275#else
276            tabres_child(:,:,:,:) = tabres(:,:,:,1:jpts)
277#endif
278         
279! WARNING :
280! tabres replaced by tabres_child in the following
281! k1 k2 replaced by 1 jpk
282! VERTICAL REFINEMENT END
283
284         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
285            ! Add asselin part
286            DO jn = n1,n2-1
287               DO jk=1,jpk
288                  DO jj=j1,j2
289                     DO ji=i1,i2
290                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
291                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
292                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
293                                 &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
294                        ENDIF
295                     ENDDO
296                  ENDDO
297               ENDDO
298            ENDDO
299         ENDIF
300         DO jn = n1,n2-1
301            DO jk=1,jpk
302               DO jj=j1,j2
303                  DO ji=i1,i2
304                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
305                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
306                     END IF
307                  END DO
308               END DO
309            END DO
310         END DO
311      ENDIF
312      !
313   END SUBROUTINE updateTS
314
315
316   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
317      !!---------------------------------------------
318      !!           *** ROUTINE updateu ***
319      !!---------------------------------------------
320      INTEGER                                 , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
321      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2), INTENT(inout) :: tabres
322      LOGICAL                                 , INTENT(in   ) :: before
323      !
324      INTEGER  ::   ji, jj, jk
325      REAL(wp) ::   zrhoy
326! VERTICAL REFINEMENT BEGIN
327      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
328      REAL(wp) :: h_in(k1:k2)
329      REAL(wp) :: h_out(1:jpk)
330      INTEGER :: N_in, N_out
331      REAL(wp) :: h_diff
332      REAL(wp) :: tabin(k1:k2)
333! VERTICAL REFINEMENT END
334      !!---------------------------------------------
335      !
336      IF( before ) THEN
337         print *, i1,i2,j1,j2,k1,k2
338         zrhoy = Agrif_Rhoy()
339# if defined key_vertical
340         AGRIF_SpecialValue = -999._wp
341         DO jk=k1,k2
342            DO jj=j1,j2
343               DO ji=i1,i2
344                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
345                                       - (umask(ji,jj,jk)-1)*999._wp
346                  tabres(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
347                                       - (umask(ji,jj,jk)-1)*999._wp
348               END DO
349            END DO
350         END DO
351#else
352            DO jk = k1,k2
353               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)
354            END DO
355#endif
356      ELSE
357         tabres_child(:,:,:) = 0.
358# if defined key_vertical
359         AGRIF_SpecialValue = -999._wp
360         DO jj=j1,j2
361         DO ji=i1,i2
362           N_in = 0
363           DO jk=k1,k2 !k2=jpk of child grid
364             IF (tabres(ji,jj,jk,2) > -900) EXIT
365             N_in = N_in + 1
366             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
367             h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
368           ENDDO
369           N_out = 0
370           DO jk=1,jpk
371             IF (umask(ji,jj,jk) == 0) EXIT
372             N_out = N_out + 1
373             h_out(N_out) = e3u_n(ji,jj,jk)
374           ENDDO
375           IF (N_in * N_out > 0) THEN
376             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
377! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
378             if (h_diff < -1.e-4) then
379             print *,'CHECK YOUR bathy U points ...',ji,jj,h_diff,e2u(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out
380             stop
381!             else ! Extends with 0
382!             N_in = N_in + 1
383!             tabin(N_in) = 0.
384!             h_in(N_in) = h_diff
385             endif
386             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)
387          ENDIF
388         ENDDO
389         ENDDO
390#else
391            DO jk=1,jpk
392               DO jj=j1,j2
393                  DO ji=i1,i2
394                     tabres_child(:,:,:) = tabres(:,:,:,1)* r1_e2u(ji,jj) / e3u_n(ji,jj,jk)
395                  END DO
396               END DO
397            END DO
398#endif
399         
400! WARNING :
401! tabres replaced by tabres_child in the following
402! k1 k2 replaced by 1 jpk
403! remove division by fs e3u (already done)
404! VERTICAL REFINEMENT END
405
406         DO jk=1,jpk
407            DO jj=j1,j2
408               DO ji=i1,i2
409!Following line now replaced by division higher up in vertical refinement case
410!                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)
411                  !
412                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
413                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
414                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
415                  ENDIF
416                  !
417                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
418               END DO
419            END DO
420         END DO
421      ENDIF
422      !
423   END SUBROUTINE updateu
424
425
426   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before )
427      !!---------------------------------------------
428      !!           *** ROUTINE updatev ***
429      !!---------------------------------------------
430      INTEGER :: i1,i2,j1,j2,k1,k2,n1,n2
431      INTEGER :: ji,jj,jk
432      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2) :: tabres
433      LOGICAL :: before
434      !!
435      REAL(wp) :: zrhox
436! VERTICAL REFINEMENT BEGIN
437      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
438      REAL(wp) :: h_in(k1:k2)
439      REAL(wp) :: h_out(1:jpk)
440      INTEGER :: N_in, N_out
441      REAL(wp) :: h_diff
442      REAL(wp) :: tabin(k1:k2)
443! VERTICAL REFINEMENT END
444      !!---------------------------------------------     
445      !
446      IF (before) THEN
447         zrhox = Agrif_Rhox()
448#if defined key_vertical
449         DO jk=k1,k2
450            DO jj=j1,j2
451               DO ji=i1,i2
452                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk)
453                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) 
454               END DO
455            END DO
456         END DO
457#else
458            DO jk=k1,k2
459               DO jj=j1,j2
460                  DO ji=i1,i2
461                     tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
462                  END DO
463               END DO
464            END DO
465#endif
466      ELSE
467         tabres_child(:,:,:) = 0.
468! VERTICAL REFINEMENT BEGIN
469#if defined key_vertical
470         DO jj=j1,j2
471         DO ji=i1,i2
472           N_in = 0
473           DO jk=k1,k2
474             IF (tabres(ji,jj,jk,2) == 0) EXIT
475             N_in = N_in + 1
476             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
477             h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
478           ENDDO
479           N_out = 0
480           DO jk=1,jpk
481             IF (vmask(ji,jj,jk) == 0) EXIT
482             N_out = N_out + 1
483             h_out(N_out) = e3v_n(ji,jj,jk)
484           ENDDO
485           IF (N_in * N_out > 0) THEN
486             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
487! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
488             if (h_diff < -1.e-4) then
489             print *,'CHECK YOUR BATHY ...'
490             stop
491!             else ! Extends with 0
492!             N_in = N_in + 1
493!             tabin(N_in) = 0.
494!             h_in(N_in) = h_diff
495             endif
496             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)
497          ENDIF
498         ENDDO
499         ENDDO
500#else
501            DO jk=k1,k2
502               DO jj=j1,j2
503                  DO ji=i1,i2
504                     tabres(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
505                  END DO
506               END DO
507            END DO
508#endif
509         
510! WARNING :
511! tabres replaced by tabres_child in the following
512! k1 k2 replaced by 1 jpk
513! remove division by fs e3v (already done)
514! VERTICAL REFINEMENT END
515
516         DO jk=k1,jpk
517            DO jj=j1,j2
518               DO ji=i1,i2
519!Following line now replaced by division higher up in vertical refinement case
520!                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)
521                  !
522                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
523                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
524                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
525                  ENDIF
526                  !
527                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
528               END DO
529            END DO
530         END DO
531      ENDIF
532      !
533   END SUBROUTINE updatev
534
535
536   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
537      !!---------------------------------------------
538      !!          *** ROUTINE updateu2d ***
539      !!---------------------------------------------
540      INTEGER, INTENT(in) :: i1, i2, j1, j2
541      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
542      LOGICAL, INTENT(in) :: before
543      !!
544      INTEGER :: ji, jj, jk
545      REAL(wp) :: zrhoy
546      REAL(wp) :: zcorr
547      !!---------------------------------------------
548      !
549      IF (before) THEN
550         zrhoy = Agrif_Rhoy()
551         DO jj=j1,j2
552            DO ji=i1,i2
553               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
554            END DO
555         END DO
556      ELSE
557         DO jj=j1,j2
558            DO ji=i1,i2
559               tabres(ji,jj) =  tabres(ji,jj) * r1_hu_n(ji,jj) * r1_e2u(ji,jj) 
560               !   
561               ! Update "now" 3d velocities:
562               spgu(ji,jj) = 0._wp
563               DO jk=1,jpkm1
564                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
565               END DO
566               spgu(ji,jj) = spgu(ji,jj) * r1_hu_n(ji,jj)
567               !
568               zcorr = tabres(ji,jj) - spgu(ji,jj)
569               DO jk=1,jpkm1             
570                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
571               END DO
572               !
573               ! Update barotropic velocities:
574               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
575                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
576                     zcorr = tabres(ji,jj) - un_b(ji,jj)
577                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
578                  END IF
579               ENDIF             
580               un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1)
581               !       
582               ! Correct "before" velocities to hold correct bt component:
583               spgu(ji,jj) = 0.e0
584               DO jk=1,jpkm1
585                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
586               END DO
587               spgu(ji,jj) = spgu(ji,jj) * r1_hu_b(ji,jj)
588               !
589               zcorr = ub_b(ji,jj) - spgu(ji,jj)
590               DO jk=1,jpkm1             
591                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
592               END DO
593               !
594            END DO
595         END DO
596      ENDIF
597      !
598   END SUBROUTINE updateu2d
599
600
601   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
602      !!---------------------------------------------
603      !!          *** ROUTINE updatev2d ***
604      !!---------------------------------------------
605      INTEGER, INTENT(in) :: i1, i2, j1, j2
606      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
607      LOGICAL, INTENT(in) :: before
608      !!
609      INTEGER :: ji, jj, jk
610      REAL(wp) :: zrhox
611      REAL(wp) :: zcorr
612      !!---------------------------------------------
613      !
614      IF (before) THEN
615         zrhox = Agrif_Rhox()
616         DO jj=j1,j2
617            DO ji=i1,i2
618               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
619            END DO
620         END DO
621      ELSE
622         DO jj=j1,j2
623            DO ji=i1,i2
624               tabres(ji,jj) =  tabres(ji,jj) * r1_hv_n(ji,jj) * r1_e1v(ji,jj) 
625               !   
626               ! Update "now" 3d velocities:
627               spgv(ji,jj) = 0.e0
628               DO jk=1,jpkm1
629                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
630               END DO
631               spgv(ji,jj) = spgv(ji,jj) * r1_hv_n(ji,jj)
632               !
633               zcorr = tabres(ji,jj) - spgv(ji,jj)
634               DO jk=1,jpkm1             
635                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
636               END DO
637               !
638               ! Update barotropic velocities:
639               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
640                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
641                     zcorr = tabres(ji,jj) - vn_b(ji,jj)
642                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
643                  END IF
644               ENDIF             
645               vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1)
646               !       
647               ! Correct "before" velocities to hold correct bt component:
648               spgv(ji,jj) = 0.e0
649               DO jk=1,jpkm1
650                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
651               END DO
652               spgv(ji,jj) = spgv(ji,jj) * r1_hv_b(ji,jj)
653               !
654               zcorr = vb_b(ji,jj) - spgv(ji,jj)
655               DO jk=1,jpkm1             
656                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
657               END DO
658               !
659            END DO
660         END DO
661      ENDIF
662      !
663   END SUBROUTINE updatev2d
664
665
666   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
667      !!---------------------------------------------
668      !!          *** ROUTINE updateSSH ***
669      !!---------------------------------------------
670      INTEGER, INTENT(in) :: i1, i2, j1, j2
671      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
672      LOGICAL, INTENT(in) :: before
673      !!
674      INTEGER :: ji, jj
675      !!---------------------------------------------
676      !
677      IF (before) THEN
678         DO jj=j1,j2
679            DO ji=i1,i2
680               tabres(ji,jj) = sshn(ji,jj)
681            END DO
682         END DO
683      ELSE
684         IF( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN
685            IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
686               DO jj=j1,j2
687                  DO ji=i1,i2
688                     sshb(ji,jj) =   sshb(ji,jj) &
689                           & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
690                  END DO
691               END DO
692            ENDIF
693         ENDIF
694         !
695         DO jj=j1,j2
696            DO ji=i1,i2
697               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
698            END DO
699         END DO
700      ENDIF
701      !
702   END SUBROUTINE updateSSH
703
704
705   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
706      !!---------------------------------------------
707      !!          *** ROUTINE updateub2b ***
708      !!---------------------------------------------
709      INTEGER, INTENT(in) :: i1, i2, j1, j2
710      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
711      LOGICAL, INTENT(in) :: before
712      !!
713      INTEGER :: ji, jj
714      REAL(wp) :: zrhoy
715      !!---------------------------------------------
716      !
717      IF (before) THEN
718         zrhoy = Agrif_Rhoy()
719         DO jj=j1,j2
720            DO ji=i1,i2
721               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
722            END DO
723         END DO
724         tabres = zrhoy * tabres
725      ELSE
726         DO jj=j1,j2
727            DO ji=i1,i2
728               ub2_b(ji,jj) = tabres(ji,jj) / e2u(ji,jj)
729            END DO
730         END DO
731      ENDIF
732      !
733   END SUBROUTINE updateub2b
734
735
736   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
737      !!---------------------------------------------
738      !!          *** ROUTINE updatevb2b ***
739      !!---------------------------------------------
740      INTEGER, INTENT(in) :: i1, i2, j1, j2
741      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
742      LOGICAL, INTENT(in) :: before
743      !!
744      INTEGER :: ji, jj
745      REAL(wp) :: zrhox
746      !!---------------------------------------------
747      !
748      IF (before) THEN
749         zrhox = Agrif_Rhox()
750         DO jj=j1,j2
751            DO ji=i1,i2
752               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
753            END DO
754         END DO
755         tabres = zrhox * tabres
756      ELSE
757         DO jj=j1,j2
758            DO ji=i1,i2
759               vb2_b(ji,jj) = tabres(ji,jj) / e1v(ji,jj)
760            END DO
761         END DO
762      ENDIF
763      !
764   END SUBROUTINE updatevb2b
765
766
767   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
768      ! currently not used
769      !!---------------------------------------------
770      !!           *** ROUTINE updateT ***
771      !!---------------------------------------------
772      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
773      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
774      LOGICAL, iNTENT(in) :: before
775      !
776      INTEGER :: ji,jj,jk
777      REAL(wp) :: ztemp
778      !!---------------------------------------------
779
780      IF (before) THEN
781         DO jk=k1,k2
782            DO jj=j1,j2
783               DO ji=i1,i2
784                  tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
785                  tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
786                  tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
787               END DO
788            END DO
789         END DO
790         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
791         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
792         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
793      ELSE
794         DO jk=k1,k2
795            DO jj=j1,j2
796               DO ji=i1,i2
797                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
798                     print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
799                     print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
800                     print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
801                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
802                     print *,'CORR = ',ztemp-1.
803                     print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
804                           tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
805                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
806                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
807                  END IF
808               END DO
809            END DO
810         END DO
811      ENDIF
812      !
813   END SUBROUTINE update_scales
814
815# if defined key_zdftke
816
817   SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
818      !!---------------------------------------------
819      !!           *** ROUTINE updateen ***
820      !!---------------------------------------------
821      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
822      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
823      LOGICAL, INTENT(in) :: before
824      !!---------------------------------------------
825      !
826      IF (before) THEN
827         ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
828      ELSE
829         en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
830      ENDIF
831      !
832   END SUBROUTINE updateEN
833
834
835   SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
836      !!---------------------------------------------
837      !!           *** ROUTINE updateavt ***
838      !!---------------------------------------------
839      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
840      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
841      LOGICAL, INTENT(in) :: before
842      !!---------------------------------------------
843      !
844      IF (before) THEN
845         ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
846      ELSE
847         avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
848      ENDIF
849      !
850   END SUBROUTINE updateAVT
851
852
853   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
854      !!---------------------------------------------
855      !!           *** ROUTINE updateavm ***
856      !!---------------------------------------------
857      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
858      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
859      LOGICAL, INTENT(in) :: before
860      !!---------------------------------------------
861      !
862      IF (before) THEN
863         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
864      ELSE
865         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
866      ENDIF
867      !
868   END SUBROUTINE updateAVM
869
870# endif /* key_zdftke */ 
871
872#else
873CONTAINS
874   SUBROUTINE agrif_opa_update_empty
875      !!---------------------------------------------
876      !!   *** ROUTINE agrif_opa_update_empty ***
877      !!---------------------------------------------
878      WRITE(*,*)  'agrif_opa_update : You should not have seen this print! error?'
879   END SUBROUTINE agrif_opa_update_empty
880#endif
881END MODULE agrif_opa_update
882
Note: See TracBrowser for help on using the repository browser.