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_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 7988

Last change on this file since 7988 was 7988, checked in by jchanut, 7 years ago

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