source: trunk/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 @ 4528

Last change on this file since 4528 was 4491, checked in by jchanut, 8 years ago

Missing Asselin correction with Agrif, #1203

  • Property svn:keywords set to Id
File size: 15.9 KB
Line 
1#define TWO_WAY
2
3MODULE agrif_opa_update
4#if defined key_agrif  && ! defined key_offline
5   USE par_oce
6   USE oce
7   USE dom_oce
8   USE agrif_oce
9   USE in_out_manager  ! I/O manager
10   USE lib_mpp
11   USE wrk_nemo 
12   USE dynspg_oce
13
14   IMPLICIT NONE
15   PRIVATE
16
17   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn
18
19   INTEGER, PUBLIC :: nbcline = 0
20
21   !!----------------------------------------------------------------------
22   !! NEMO/NST 3.3 , NEMO Consortium (2010)
23   !! $Id$
24   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
25   !!----------------------------------------------------------------------
26
27CONTAINS
28
29   SUBROUTINE Agrif_Update_Tra( kt )
30      !!---------------------------------------------
31      !!   *** ROUTINE Agrif_Update_Tra ***
32      !!---------------------------------------------
33      !!
34      INTEGER, INTENT(in) :: kt
35      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztab
36
37
38      IF((Agrif_NbStepint() .NE. (Agrif_irhot()-1)).AND.(kt /= 0)) RETURN
39#if defined TWO_WAY
40      CALL wrk_alloc( jpi, jpj, jpk, jpts, ztab )
41
42      Agrif_UseSpecialValueInUpdate = .TRUE.
43      Agrif_SpecialValueFineGrid = 0.
44
45      IF (MOD(nbcline,nbclineupdate) == 0) THEN
46         CALL Agrif_Update_Variable(ztab,tsn_id, procname=updateTS)
47      ELSE
48         CALL Agrif_Update_Variable(ztab,tsn_id,locupdate=(/0,2/), procname=updateTS)
49      ENDIF
50
51      Agrif_UseSpecialValueInUpdate = .FALSE.
52
53      CALL wrk_dealloc( jpi, jpj, jpk, jpts, ztab )
54#endif
55
56   END SUBROUTINE Agrif_Update_Tra
57
58   SUBROUTINE Agrif_Update_Dyn( kt )
59      !!---------------------------------------------
60      !!   *** ROUTINE Agrif_Update_Dyn ***
61      !!---------------------------------------------
62      !!
63      INTEGER, INTENT(in) :: kt
64      REAL(wp), POINTER, DIMENSION(:,:) :: ztab2d
65      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztab
66
67
68      IF ((Agrif_NbStepint() .NE. (Agrif_irhot()-1)).AND.(kt /= 0)) Return
69#if defined TWO_WAY
70      CALL wrk_alloc( jpi, jpj,      ztab2d )
71      CALL wrk_alloc( jpi, jpj, jpk, ztab   )
72
73      IF (mod(nbcline,nbclineupdate) == 0) THEN
74         CALL Agrif_Update_Variable(ztab,un_id,procname = updateU)
75         CALL Agrif_Update_Variable(ztab,vn_id,procname = updateV)
76      ELSE
77         CALL Agrif_Update_Variable(ztab,un_id,locupdate=(/0,1/),procname = updateU)
78         CALL Agrif_Update_Variable(ztab,vn_id,locupdate=(/0,1/),procname = updateV)         
79      ENDIF
80
81      CALL Agrif_Update_Variable(ztab2d,e1u_id,procname = updateU2d)
82      CALL Agrif_Update_Variable(ztab2d,e2v_id,procname = updateV2d)
83
84#if defined key_dynspg_ts
85      IF (ln_bt_fw) THEN
86         ! Update time integrated transports
87         IF (mod(nbcline,nbclineupdate) == 0) THEN
88            CALL Agrif_Update_Variable(ztab2d,ub2b_id,procname = updateub2b)
89            CALL Agrif_Update_Variable(ztab2d,vb2b_id,procname = updatevb2b)
90         ELSE
91            CALL Agrif_Update_Variable(ztab2d,ub2b_id,locupdate=(/0,1/),procname = updateub2b)
92            CALL Agrif_Update_Variable(ztab2d,vb2b_id,locupdate=(/0,1/),procname = updatevb2b)
93         ENDIF
94      END IF 
95#endif
96
97      nbcline = nbcline + 1
98
99      Agrif_UseSpecialValueInUpdate = .TRUE. 
100      Agrif_SpecialValueFineGrid = 0.
101      CALL Agrif_Update_Variable(ztab2d,sshn_id,procname = updateSSH)
102      Agrif_UseSpecialValueInUpdate = .FALSE.
103
104      CALL wrk_dealloc( jpi, jpj,      ztab2d )
105      CALL wrk_dealloc( jpi, jpj, jpk, ztab   )
106
107!Done in step
108!      CALL Agrif_ChildGrid_To_ParentGrid()
109!      CALL recompute_diags( kt )
110!      CALL Agrif_ParentGrid_To_ChildGrid()
111
112#endif
113
114   END SUBROUTINE Agrif_Update_Dyn
115
116   SUBROUTINE recompute_diags( kt )
117      !!---------------------------------------------
118      !!   *** ROUTINE recompute_diags ***
119      !!---------------------------------------------
120      INTEGER, INTENT(in) :: kt
121
122   END SUBROUTINE recompute_diags
123
124   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
125      !!---------------------------------------------
126      !!           *** ROUTINE updateT ***
127      !!---------------------------------------------
128#  include "domzgr_substitute.h90"
129
130      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
131      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
132      LOGICAL, iNTENT(in) :: before
133
134      INTEGER :: ji,jj,jk,jn
135
136      IF (before) THEN
137         DO jn = n1,n2
138            DO jk=k1,k2
139               DO jj=j1,j2
140                  DO ji=i1,i2
141                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
142                  END DO
143               END DO
144            END DO
145         END DO
146      ELSE
147         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
148         ! Add asselin part
149            DO jn = n1,n2
150               DO jk=k1,k2
151                  DO jj=j1,j2
152                     DO ji=i1,i2
153                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
154                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
155                              & + atfp * ( tabres(ji,jj,jk,jn) &
156                              &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
157                        ENDIF
158                     ENDDO
159                  ENDDO
160               ENDDO
161            ENDDO
162         ENDIF
163
164         DO jn = n1,n2
165            DO jk=k1,k2
166               DO jj=j1,j2
167                  DO ji=i1,i2
168                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
169                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk)
170                     END IF
171                  END DO
172               END DO
173            END DO
174         END DO
175      ENDIF
176
177   END SUBROUTINE updateTS
178
179   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before )
180      !!---------------------------------------------
181      !!           *** ROUTINE updateu ***
182      !!---------------------------------------------
183#  include "domzgr_substitute.h90"
184
185      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
186      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
187      LOGICAL, INTENT(in) :: before
188
189      INTEGER :: ji, jj, jk
190      REAL(wp) :: zrhoy
191
192      IF (before) THEN
193         zrhoy = Agrif_Rhoy()
194         DO jk=k1,k2
195            DO jj=j1,j2
196               DO ji=i1,i2
197                  tabres(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk)
198                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3u_n(ji,jj,jk)
199               END DO
200            END DO
201         END DO
202         tabres = zrhoy * tabres
203      ELSE
204         DO jk=k1,k2
205            DO jj=j1,j2
206               DO ji=i1,i2
207                  tabres(ji,jj,jk) = tabres(ji,jj,jk) / e2u(ji,jj) / fse3u_n(ji,jj,jk)
208                  !
209                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
210                     ub(ji,jj,jk) = ub(ji,jj,jk) & 
211                       & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
212                  ENDIF
213                  !
214                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk)
215               END DO
216            END DO
217         END DO
218      ENDIF
219
220   END SUBROUTINE updateu
221
222   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before )
223      !!---------------------------------------------
224      !!           *** ROUTINE updatev ***
225      !!---------------------------------------------
226#  include "domzgr_substitute.h90"
227
228      INTEGER :: i1,i2,j1,j2,k1,k2
229      INTEGER :: ji,jj,jk
230      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: tabres
231      LOGICAL :: before
232
233      REAL(wp) :: zrhox
234
235      IF (before) THEN
236         zrhox = Agrif_Rhox()
237         DO jk=k1,k2
238            DO jj=j1,j2
239               DO ji=i1,i2
240                  tabres(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk)
241                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3v_n(ji,jj,jk)
242               END DO
243            END DO
244         END DO
245         tabres = zrhox * tabres
246      ELSE
247         DO jk=k1,k2
248            DO jj=j1,j2
249               DO ji=i1,i2
250                  tabres(ji,jj,jk) = tabres(ji,jj,jk) / e1v(ji,jj) / fse3v_n(ji,jj,jk)
251                  !
252                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
253                     vb(ji,jj,jk) = vb(ji,jj,jk) & 
254                       & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
255                  ENDIF
256                  !
257                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk)
258               END DO
259            END DO
260         END DO
261      ENDIF
262
263   END SUBROUTINE updatev
264
265   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
266      !!---------------------------------------------
267      !!          *** ROUTINE updateu2d ***
268      !!---------------------------------------------
269#  include "domzgr_substitute.h90"
270
271      INTEGER, INTENT(in) :: i1, i2, j1, j2
272      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
273      LOGICAL, INTENT(in) :: before
274
275      INTEGER :: ji, jj, jk
276      REAL(wp) :: zrhoy
277      REAL(wp) :: zcorr
278
279      IF (before) THEN
280         zrhoy = Agrif_Rhoy()
281         DO jj=j1,j2
282            DO ji=i1,i2
283               tabres(ji,jj) = un_b(ji,jj) * hu(ji,jj) * e2u(ji,jj)
284            END DO
285         END DO
286         tabres = zrhoy * tabres
287      ELSE
288         DO jj=j1,j2
289            DO ji=i1,i2
290               tabres(ji,jj) =  tabres(ji,jj) * hur(ji,jj) / e2u(ji,jj) 
291               !   
292               ! Update "now" 3d velocities:
293               spgu(ji,jj) = 0.e0
294               DO jk=1,jpkm1
295                  spgu(ji,jj) = spgu(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk)
296               END DO
297               spgu(ji,jj) = spgu(ji,jj) * hur(ji,jj)
298               !
299               zcorr = tabres(ji,jj) - spgu(ji,jj)
300               DO jk=1,jpkm1             
301                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
302               END DO
303               !
304               ! Update barotropic velocities:
305#if defined key_dynspg_ts
306               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
307                  zcorr = tabres(ji,jj) - un_b(ji,jj)
308                  ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1)
309               END IF
310#endif               
311               un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1)
312               !       
313               ! Correct "before" velocities to hold correct bt component:
314               spgu(ji,jj) = 0.e0
315               DO jk=1,jpkm1
316                  spgu(ji,jj) = spgu(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)
317               END DO
318               spgu(ji,jj) = spgu(ji,jj) * hur_b(ji,jj)
319               !
320               zcorr = ub_b(ji,jj) - spgu(ji,jj)
321               DO jk=1,jpkm1             
322                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)           
323               END DO
324               !
325            END DO
326         END DO
327      ENDIF
328
329   END SUBROUTINE updateu2d
330
331   SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
332      !!---------------------------------------------
333      !!          *** ROUTINE updatev2d ***
334      !!---------------------------------------------
335
336      INTEGER, INTENT(in) :: i1, i2, j1, j2
337      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
338      LOGICAL, INTENT(in) :: before
339
340      INTEGER :: ji, jj, jk
341      REAL(wp) :: zrhox
342      REAL(wp) :: zcorr
343
344      IF (before) THEN
345         zrhox = Agrif_Rhox()
346         DO jj=j1,j2
347            DO ji=i1,i2
348               tabres(ji,jj) = vn_b(ji,jj) * hv(ji,jj) * e1v(ji,jj) 
349            END DO
350         END DO
351         tabres = zrhox * tabres
352      ELSE
353         DO jj=j1,j2
354            DO ji=i1,i2
355               tabres(ji,jj) =  tabres(ji,jj) * hvr(ji,jj) / e1v(ji,jj) 
356               !   
357               ! Update "now" 3d velocities:
358               spgv(ji,jj) = 0.e0
359               DO jk=1,jpkm1
360                  spgv(ji,jj) = spgv(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk)
361               END DO
362               spgv(ji,jj) = spgv(ji,jj) * hvr(ji,jj)
363               !
364               zcorr = tabres(ji,jj) - spgv(ji,jj)
365               DO jk=1,jpkm1             
366                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
367               END DO
368               !
369               ! Update barotropic velocities:
370#if defined key_dynspg_ts
371               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
372                  zcorr = tabres(ji,jj) - vn_b(ji,jj)
373                  vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1)
374               END IF
375#endif               
376               vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1)
377               !       
378               ! Correct "before" velocities to hold correct bt component:
379               spgv(ji,jj) = 0.e0
380               DO jk=1,jpkm1
381                  spgv(ji,jj) = spgv(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)
382               END DO
383               spgv(ji,jj) = spgv(ji,jj) * hvr_b(ji,jj)
384               !
385               zcorr = vb_b(ji,jj) - spgv(ji,jj)
386               DO jk=1,jpkm1             
387                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)           
388               END DO
389               !
390            END DO
391         END DO
392      ENDIF
393
394   END SUBROUTINE updatev2d
395
396   SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
397      !!---------------------------------------------
398      !!          *** ROUTINE updateSSH ***
399      !!---------------------------------------------
400#  include "domzgr_substitute.h90"
401
402      INTEGER, INTENT(in) :: i1, i2, j1, j2
403      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
404      LOGICAL, INTENT(in) :: before
405
406      INTEGER :: ji, jj
407
408      IF (before) THEN
409         DO jj=j1,j2
410            DO ji=i1,i2
411               tabres(ji,jj) = sshn(ji,jj)
412            END DO
413         END DO
414      ELSE
415
416#if ! defined key_dynspg_ts
417         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
418            DO jj=j1,j2
419               DO ji=i1,i2
420                sshb(ji,jj) =   sshb(ji,jj) &
421                 & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
422               END DO
423            END DO
424         ENDIF
425#endif
426         DO jj=j1,j2
427            DO ji=i1,i2
428               sshn(ji,jj) = tabres(ji,jj) * tmask(ji,jj,1)
429            END DO
430         END DO
431      ENDIF
432
433   END SUBROUTINE updateSSH
434
435   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
436      !!---------------------------------------------
437      !!          *** ROUTINE updateub2b ***
438      !!---------------------------------------------
439
440      INTEGER, INTENT(in) :: i1, i2, j1, j2
441      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
442      LOGICAL, INTENT(in) :: before
443
444      INTEGER :: ji, jj
445      REAL(wp) :: zrhoy
446
447      IF (before) THEN
448         zrhoy = Agrif_Rhoy()
449         DO jj=j1,j2
450            DO ji=i1,i2
451               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
452            END DO
453         END DO
454         tabres = zrhoy * tabres
455      ELSE
456         DO jj=j1,j2
457            DO ji=i1,i2
458               ub2_b(ji,jj) = tabres(ji,jj) / e2u(ji,jj)
459            END DO
460         END DO
461      ENDIF
462
463   END SUBROUTINE updateub2b
464
465   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
466      !!---------------------------------------------
467      !!          *** ROUTINE updatevb2b ***
468      !!---------------------------------------------
469
470      INTEGER, INTENT(in) :: i1, i2, j1, j2
471      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
472      LOGICAL, INTENT(in) :: before
473
474      INTEGER :: ji, jj
475      REAL(wp) :: zrhox
476
477      IF (before) THEN
478         zrhox = Agrif_Rhox()
479         DO jj=j1,j2
480            DO ji=i1,i2
481               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
482            END DO
483         END DO
484         tabres = zrhox * tabres
485      ELSE
486         DO jj=j1,j2
487            DO ji=i1,i2
488               vb2_b(ji,jj) = tabres(ji,jj) / e1v(ji,jj)
489            END DO
490         END DO
491      ENDIF
492
493   END SUBROUTINE updatevb2b
494
495#else
496CONTAINS
497   SUBROUTINE agrif_opa_update_empty
498      !!---------------------------------------------
499      !!   *** ROUTINE agrif_opa_update_empty ***
500      !!---------------------------------------------
501      WRITE(*,*)  'agrif_opa_update : You should not have seen this print! error?'
502   END SUBROUTINE agrif_opa_update_empty
503#endif
504END MODULE agrif_opa_update
505
Note: See TracBrowser for help on using the repository browser.