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/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 6401

Last change on this file since 6401 was 6401, checked in by timgraham, 8 years ago

Reinstate svn keywords before update to head of trunk

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