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 @ 8902

Last change on this file since 8902 was 8902, checked in by timgraham, 6 years ago

Bug fixes for non key_vertical case, reverted to interpolating tracer rather than tracer content
  • Property svn:keywords set to Id
File size: 31.0 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      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
190      REAL(wp) :: h_in(k1:k2)
191      REAL(wp) :: h_out(1:jpk)
192      INTEGER  :: N_in, N_out
193      REAL(wp) :: h_diff
194      REAL(wp) :: zrho_xy
195      REAL(wp) :: tabin(k1:k2,n1:n2)
196      !!---------------------------------------------
197      !
198      IF (before) THEN
199# if defined key_vertical
200         AGRIF_SpecialValue = -999._wp
201         zrho_xy = Agrif_rhox() * Agrif_rhoy() 
202         DO jn = n1,n2-1
203            DO jk=k1,k2
204               DO jj=j1,j2
205                  DO ji=i1,i2
206                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) &
207                                           * tmask(ji,jj,jk) * zrho_xy + (tmask(ji,jj,jk)-1)*999._wp
208                  END DO
209               END DO
210            END DO
211         END DO
212         DO jk=k1,k2
213            DO jj=j1,j2
214               DO ji=i1,i2
215                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy  &
216                                           + (tmask(ji,jj,jk)-1)*999._wp
217               END DO
218            END DO
219         END DO
220#else
221         DO jn = n1,n2-1
222            DO jk=k1,k2
223               DO jj=j1,j2
224                  DO ji=i1,i2
225                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
226                  END DO
227               END DO
228            END DO
229         END DO       
230#endif
231      ELSE
232         tabres_child(:,:,:,:) = 0.
233# if defined key_vertical
234         AGRIF_SpecialValue = 0._wp
235         DO jj=j1,j2
236            DO ji=i1,i2
237               N_in = 0
238               DO jk=k1,k2 !k2 = jpk of child grid
239                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT
240                  N_in = N_in + 1
241                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
242                  h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj)
243               ENDDO
244               N_out = 0
245               DO jk=1,jpk ! jpk of parent grid
246                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
247                  N_out = N_out + 1
248                  h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
249               ENDDO
250               IF (N_in > 0) THEN !Remove this?
251                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
252                  IF (h_diff < -1.e-4) THEN
253                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
254                     print *,h_in(1:N_in)
255                     print *,h_out(1:N_out)
256                     STOP
257                  ENDIF
258                  DO jn=n1,n2-1
259                     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)
260                  ENDDO
261               ENDIF
262            ENDDO
263         ENDDO
264#else
265         tabres_child(:,:,:,:) = tabres(:,:,:,1:jpts)
266#endif
267         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
268            ! Add asselin part
269            DO jn = n1,n2-1
270               DO jk=1,jpk
271                  DO jj=j1,j2
272                     DO ji=i1,i2
273                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
274                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
275                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
276                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
277                        ENDIF
278                     ENDDO
279                  ENDDO
280               ENDDO
281            ENDDO
282         ENDIF
283         DO jn = n1,n2-1
284            DO jk=1,jpk
285               DO jj=j1,j2
286                  DO ji=i1,i2
287                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
288                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
289                     END IF
290                  END DO
291               END DO
292            END DO
293         END DO
294      ENDIF
295      !
296   END SUBROUTINE updateTS
297
298
299   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
300      !!---------------------------------------------
301      !!           *** ROUTINE updateu ***
302      !!---------------------------------------------
303      INTEGER                                 , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
304      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2), INTENT(inout) :: tabres
305      LOGICAL                                 , INTENT(in   ) :: before
306      !
307      INTEGER  ::   ji, jj, jk
308      REAL(wp) ::   zrhoy
309! VERTICAL REFINEMENT BEGIN
310      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
311      REAL(wp) :: h_in(k1:k2)
312      REAL(wp) :: h_out(1:jpk)
313      INTEGER :: N_in, N_out
314      REAL(wp) :: h_diff, excess, thick
315      REAL(wp) :: tabin(k1:k2)
316! VERTICAL REFINEMENT END
317      !!---------------------------------------------
318      !
319      IF( before ) THEN
320         zrhoy = Agrif_Rhoy()
321# if defined key_vertical
322         AGRIF_SpecialValue = -999._wp
323         DO jk=k1,k2
324            DO jj=j1,j2
325               DO ji=i1,i2
326                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
327                                       + (umask(ji,jj,jk)-1)*999._wp
328                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
329                                       + (umask(ji,jj,jk)-1)*999._wp
330               END DO
331            END DO
332         END DO
333#else
334         DO jk = k1,k2
335            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)
336         END DO
337#endif
338      ELSE
339         tabres_child(:,:,:) = 0.
340# if defined key_vertical
341         AGRIF_SpecialValue = 0._wp
342         DO jj=j1,j2
343            DO ji=i1,i2
344               N_in = 0
345               h_in(:) = 0._wp
346               tabin(:) = 0._wp
347               DO jk=k1,k2 !k2=jpk of child grid
348                  IF( tabres(ji,jj,jk,2) < -900) EXIT
349                  N_in = N_in + 1
350                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
351                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
352               ENDDO
353               N_out = 0
354               DO jk=1,jpk
355                  IF (umask(ji,jj,jk) == 0) EXIT
356                  N_out = N_out + 1
357                  h_out(N_out) = e3u_n(ji,jj,jk)
358               ENDDO
359               IF (N_in * N_out > 0) THEN
360                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
361                  IF (h_diff < -1.e-4) THEN
362!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
363!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
364                     excess = 0._wp
365                     DO jk=N_in,1,-1
366                        thick = MIN(-1*h_diff, h_in(jk))
367                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
368                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
369                        h_diff = h_diff + thick
370                        IF ( h_diff == 0) THEN
371                           N_in = jk
372                           h_in(jk) = h_in(jk) - thick
373                           EXIT
374                        ENDIF
375                     ENDDO
376                  ENDIF
377                  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)
378                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
379               ENDIF
380            ENDDO
381         ENDDO
382#else
383            DO jk=1,jpk
384               DO jj=j1,j2
385                  DO ji=i1,i2
386                     tabres_child(:,:,:) = tabres(:,:,:,1)* r1_e2u(ji,jj) / e3u_n(ji,jj,jk)
387                  END DO
388               END DO
389            END DO
390#endif
391
392         DO jk=1,jpk
393            DO jj=j1,j2
394               DO ji=i1,i2
395                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
396                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
397                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
398                  ENDIF
399                  !
400                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
401               END DO
402            END DO
403         END DO
404      ENDIF
405      !
406   END SUBROUTINE updateu
407
408
409   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before )
410      !!---------------------------------------------
411      !!           *** ROUTINE updatev ***
412      !!---------------------------------------------
413      INTEGER :: i1,i2,j1,j2,k1,k2,n1,n2
414      INTEGER :: ji,jj,jk
415      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2) :: tabres
416      LOGICAL :: before
417      !!
418      REAL(wp) :: zrhox
419! VERTICAL REFINEMENT BEGIN
420      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
421      REAL(wp) :: h_in(k1:k2)
422      REAL(wp) :: h_out(1:jpk)
423      INTEGER :: N_in, N_out
424      REAL(wp) :: h_diff, excess, thick
425      REAL(wp) :: tabin(k1:k2)
426! VERTICAL REFINEMENT END
427      !!---------------------------------------------     
428      !
429      IF (before) THEN
430         zrhox = Agrif_Rhox()
431#if defined key_vertical
432         AGRIF_SpecialValue = -999._wp
433         DO jk=k1,k2
434            DO jj=j1,j2
435               DO ji=i1,i2
436                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
437                                       + (vmask(ji,jj,jk)-1)*999._wp
438                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
439                                       + (vmask(ji,jj,jk)-1)*999._wp
440               END DO
441            END DO
442         END DO
443#else
444            DO jk=k1,k2
445               DO jj=j1,j2
446                  DO ji=i1,i2
447                     tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
448                  END DO
449               END DO
450            END DO
451#endif
452      ELSE
453         tabres_child(:,:,:) = 0.
454#if defined key_vertical
455         AGRIF_SpecialValue = 0._wp
456         DO jj=j1,j2
457            DO ji=i1,i2
458               N_in = 0
459               DO jk=k1,k2
460                  IF (tabres(ji,jj,jk,2) < -900) EXIT
461                  N_in = N_in + 1
462                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
463                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
464               ENDDO
465               N_out = 0
466               DO jk=1,jpk
467                  IF (vmask(ji,jj,jk) == 0) EXIT
468                  N_out = N_out + 1
469                  h_out(N_out) = e3v_n(ji,jj,jk)
470               ENDDO
471               IF (N_in * N_out > 0) THEN
472                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
473                  IF (h_diff < -1.e-4) then
474!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.
475!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
476                     excess = 0._wp
477                     DO jk=N_in,1,-1
478                        thick = MIN(-1*h_diff, h_in(jk))
479                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
480                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
481                        h_diff = h_diff + thick
482                        IF ( h_diff == 0) THEN
483                           N_in = jk
484                           h_in(jk) = h_in(jk) - thick
485                           EXIT
486                        ENDIF
487                     ENDDO
488                  ENDIF
489                  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)
490                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
491               ENDIF
492            ENDDO
493         ENDDO
494#else
495         DO jk=k1,k2
496            DO jj=j1,j2
497               DO ji=i1,i2
498                  tabres(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
499               END DO
500            END DO
501         END DO
502#endif
503
504         DO jk=k1,jpk
505            DO jj=j1,j2
506               DO ji=i1,i2
507                  !
508                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
509                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
510                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
511                  ENDIF
512                  !
513                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
514               END DO
515            END DO
516         END DO
517      ENDIF
518      !
519   END SUBROUTINE updatev
520
521
522   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
523      !!---------------------------------------------
524      !!          *** ROUTINE updateu2d ***
525      !!---------------------------------------------
526      INTEGER, INTENT(in) :: i1, i2, j1, j2
527      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
528      LOGICAL, INTENT(in) :: before
529      !!
530      INTEGER :: ji, jj, jk
531      REAL(wp) :: zrhoy
532      REAL(wp) :: zcorr
533      !!---------------------------------------------
534      !
535      IF (before) THEN
536         zrhoy = Agrif_Rhoy()
537         DO jj=j1,j2
538            DO ji=i1,i2
539               tabres(ji,jj) = zrhoy * un_b(ji,jj) * hu_n(ji,jj) * e2u(ji,jj)
540            END DO
541         END DO
542      ELSE
543         DO jj=j1,j2
544            DO ji=i1,i2
545               tabres(ji,jj) =  tabres(ji,jj) * r1_hu_n(ji,jj) * r1_e2u(ji,jj) 
546               !   
547               ! Update "now" 3d velocities:
548               spgu(ji,jj) = 0._wp
549               DO jk=1,jpkm1
550                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk)
551               END DO
552               spgu(ji,jj) = spgu(ji,jj) * r1_hu_n(ji,jj)
553               !
554               zcorr = tabres(ji,jj) - spgu(ji,jj)
555               DO jk=1,jpkm1             
556                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
557               END DO
558               !
559               ! Update barotropic velocities:
560               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
561                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
562                     zcorr = tabres(ji,jj) - un_b(ji,jj)
563                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
564                  END IF
565               ENDIF             
566               un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1)
567               !       
568               ! Correct "before" velocities to hold correct bt component:
569               spgu(ji,jj) = 0.e0
570               DO jk=1,jpkm1
571                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk)
572               END DO
573               spgu(ji,jj) = spgu(ji,jj) * r1_hu_b(ji,jj)
574               !
575               zcorr = ub_b(ji,jj) - spgu(ji,jj)
576               DO jk=1,jpkm1             
577                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
578               END DO
579               !
580            END DO
581         END DO
582      ENDIF
583      !
584   END SUBROUTINE updateu2d
585
586
587   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
588      !!---------------------------------------------
589      !!          *** ROUTINE updatev2d ***
590      !!---------------------------------------------
591      INTEGER, INTENT(in) :: i1, i2, j1, j2
592      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
593      LOGICAL, INTENT(in) :: before
594      !!
595      INTEGER :: ji, jj, jk
596      REAL(wp) :: zrhox
597      REAL(wp) :: zcorr
598      !!---------------------------------------------
599      !
600      IF (before) THEN
601         zrhox = Agrif_Rhox()
602         DO jj=j1,j2
603            DO ji=i1,i2
604               tabres(ji,jj) = zrhox * vn_b(ji,jj) * hv_n(ji,jj) * e1v(ji,jj) 
605            END DO
606         END DO
607      ELSE
608         DO jj=j1,j2
609            DO ji=i1,i2
610               tabres(ji,jj) =  tabres(ji,jj) * r1_hv_n(ji,jj) * r1_e1v(ji,jj) 
611               !   
612               ! Update "now" 3d velocities:
613               spgv(ji,jj) = 0.e0
614               DO jk=1,jpkm1
615                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk)
616               END DO
617               spgv(ji,jj) = spgv(ji,jj) * r1_hv_n(ji,jj)
618               !
619               zcorr = tabres(ji,jj) - spgv(ji,jj)
620               DO jk=1,jpkm1             
621                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
622               END DO
623               !
624               ! Update barotropic velocities:
625               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
626                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
627                     zcorr = tabres(ji,jj) - vn_b(ji,jj)
628                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
629                  END IF
630               ENDIF             
631               vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1)
632               !       
633               ! Correct "before" velocities to hold correct bt component:
634               spgv(ji,jj) = 0.e0
635               DO jk=1,jpkm1
636                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk)
637               END DO
638               spgv(ji,jj) = spgv(ji,jj) * r1_hv_b(ji,jj)
639               !
640               zcorr = vb_b(ji,jj) - spgv(ji,jj)
641               DO jk=1,jpkm1             
642                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
643               END DO
644               !
645            END DO
646         END DO
647      ENDIF
648      !
649   END SUBROUTINE updatev2d
650
651
652   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
653      !!---------------------------------------------
654      !!          *** ROUTINE updateSSH ***
655      !!---------------------------------------------
656      INTEGER, INTENT(in) :: i1, i2, j1, j2
657      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
658      LOGICAL, INTENT(in) :: before
659      !!
660      INTEGER :: ji, jj
661      !!---------------------------------------------
662      !
663      IF (before) THEN
664         DO jj=j1,j2
665            DO ji=i1,i2
666               tabres(ji,jj) = sshn(ji,jj)
667            END DO
668         END DO
669      ELSE
670         IF( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN
671            IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
672               DO jj=j1,j2
673                  DO ji=i1,i2
674                     sshb(ji,jj) =   sshb(ji,jj) &
675                           & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
676                  END DO
677               END DO
678            ENDIF
679         ENDIF
680         !
681         DO jj=j1,j2
682            DO ji=i1,i2
683               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
684            END DO
685         END DO
686      ENDIF
687      !
688   END SUBROUTINE updateSSH
689
690
691   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
692      !!---------------------------------------------
693      !!          *** ROUTINE updateub2b ***
694      !!---------------------------------------------
695      INTEGER, INTENT(in) :: i1, i2, j1, j2
696      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
697      LOGICAL, INTENT(in) :: before
698      !!
699      INTEGER :: ji, jj
700      REAL(wp) :: zrhoy
701      !!---------------------------------------------
702      !
703      IF (before) THEN
704         zrhoy = Agrif_Rhoy()
705         DO jj=j1,j2
706            DO ji=i1,i2
707               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
708            END DO
709         END DO
710         tabres = zrhoy * tabres
711      ELSE
712         DO jj=j1,j2
713            DO ji=i1,i2
714               ub2_b(ji,jj) = tabres(ji,jj) / e2u(ji,jj)
715            END DO
716         END DO
717      ENDIF
718      !
719   END SUBROUTINE updateub2b
720
721
722   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
723      !!---------------------------------------------
724      !!          *** ROUTINE updatevb2b ***
725      !!---------------------------------------------
726      INTEGER, INTENT(in) :: i1, i2, j1, j2
727      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
728      LOGICAL, INTENT(in) :: before
729      !!
730      INTEGER :: ji, jj
731      REAL(wp) :: zrhox
732      !!---------------------------------------------
733      !
734      IF (before) THEN
735         zrhox = Agrif_Rhox()
736         DO jj=j1,j2
737            DO ji=i1,i2
738               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
739            END DO
740         END DO
741         tabres = zrhox * tabres
742      ELSE
743         DO jj=j1,j2
744            DO ji=i1,i2
745               vb2_b(ji,jj) = tabres(ji,jj) / e1v(ji,jj)
746            END DO
747         END DO
748      ENDIF
749   END SUBROUTINE updatevb2b
750
751
752   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
753      ! currently not used
754      !!---------------------------------------------
755      !!           *** ROUTINE updateT ***
756      !!---------------------------------------------
757      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
758      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
759      LOGICAL, iNTENT(in) :: before
760      !
761      INTEGER :: ji,jj,jk
762      REAL(wp) :: ztemp
763      !!---------------------------------------------
764
765      IF (before) THEN
766         DO jk=k1,k2
767            DO jj=j1,j2
768               DO ji=i1,i2
769                  tabres(ji,jj,jk,1) = e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
770                  tabres(ji,jj,jk,2) = e1t(ji,jj)*tmask(ji,jj,jk)
771                  tabres(ji,jj,jk,3) = e2t(ji,jj)*tmask(ji,jj,jk)
772               END DO
773            END DO
774         END DO
775         tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()
776         tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox()
777         tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy()
778      ELSE
779         DO jk=k1,k2
780            DO jj=j1,j2
781               DO ji=i1,i2
782                  IF( tabres(ji,jj,jk,1) .NE. 0. ) THEN
783                     print *,'VAL = ',ji,jj,jk,tabres(ji,jj,jk,1),e1t(ji,jj)*e2t(ji,jj)*tmask(ji,jj,jk)
784                     print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk)
785                     print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk)
786                     ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))
787                     print *,'CORR = ',ztemp-1.
788                     print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, &
789                           tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp
790                     e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp
791                     e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp
792                  END IF
793               END DO
794            END DO
795         END DO
796      ENDIF
797      !
798   END SUBROUTINE update_scales
799
800# if defined key_zdftke
801
802   SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
803      !!---------------------------------------------
804      !!           *** ROUTINE updateen ***
805      !!---------------------------------------------
806      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
807      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
808      LOGICAL, INTENT(in) :: before
809      !!---------------------------------------------
810      !
811      IF (before) THEN
812         ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
813      ELSE
814         en(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
815      ENDIF
816      !
817   END SUBROUTINE updateEN
818
819
820   SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
821      !!---------------------------------------------
822      !!           *** ROUTINE updateavt ***
823      !!---------------------------------------------
824      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
825      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
826      LOGICAL, INTENT(in) :: before
827      !!---------------------------------------------
828      !
829      IF (before) THEN
830         ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
831      ELSE
832         avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
833      ENDIF
834      !
835   END SUBROUTINE updateAVT
836
837
838   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before )
839      !!---------------------------------------------
840      !!           *** ROUTINE updateavm ***
841      !!---------------------------------------------
842      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
843      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
844      LOGICAL, INTENT(in) :: before
845      !!---------------------------------------------
846      !
847      IF (before) THEN
848         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
849      ELSE
850         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
851      ENDIF
852      !
853   END SUBROUTINE updateAVM
854
855# endif /* key_zdftke */ 
856
857#else
858CONTAINS
859   SUBROUTINE agrif_opa_update_empty
860      !!---------------------------------------------
861      !!   *** ROUTINE agrif_opa_update_empty ***
862      !!---------------------------------------------
863      WRITE(*,*)  'agrif_opa_update : You should not have seen this print! error?'
864   END SUBROUTINE agrif_opa_update_empty
865#endif
866END MODULE agrif_opa_update
867
Note: See TracBrowser for help on using the repository browser.