source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_top_update.F90 @ 11463

Last change on this file since 11463 was 11463, checked in by acc, 14 months ago

Branch: dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap. Minor bugfix in step.F90 to enable AGRIF SETTE tests to run. Also merged prettification changes to NST routines from the dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch. Neither AGRIF_DEMO nor VORTEX restart perfectly (drifting after 8 and 121 timesteps, respectively).

  • Property svn:keywords set to Id
File size: 9.3 KB
Line 
1#define TWO_WAY
2#undef DECAL_FEEDBACK
3
4MODULE agrif_top_update
5   !!======================================================================
6   !!                ***  MODULE agrif_top_update  ***
7   !! AGRIF :   update package for passive tracers (TOP)
8   !!======================================================================
9   !! History : 
10   !!----------------------------------------------------------------------
11#if defined key_agrif && defined key_top
12   !!----------------------------------------------------------------------
13   !!   'key_agrif'                                              AGRIF zoom
14   !!   'key_TOP'                                           on-line tracers
15   !!----------------------------------------------------------------------
16   USE par_oce
17   USE oce
18   USE dom_oce
19   USE agrif_oce
20   USE par_trc
21   USE trc
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC Agrif_Update_Trc
27
28   !!----------------------------------------------------------------------
29   !! NEMO/NST 4.0 , NEMO Consortium (2018)
30   !! $Id$
31   !! Software governed by the CeCILL license (see ./LICENSE)
32   !!----------------------------------------------------------------------
33CONTAINS
34
35   SUBROUTINE Agrif_Update_Trc( )
36      !!----------------------------------------------------------------------
37      !!                   *** ROUTINE Agrif_Update_Trc ***
38      !!----------------------------------------------------------------------
39      !
40      IF (Agrif_Root()) RETURN 
41      !
42#if defined TWO_WAY   
43      Agrif_UseSpecialValueInUpdate = .TRUE.
44      Agrif_SpecialValueFineGrid    = 0._wp
45      !
46# if ! defined DECAL_FEEDBACK
47      CALL Agrif_Update_Variable(trn_id, procname=updateTRC )
48!      CALL Agrif_Update_Variable( trn_id, locupdate=(/0,2/), procname=updateTRC )
49# else
50      CALL Agrif_Update_Variable(trn_id, locupdate=(/1,0/),procname=updateTRC )
51!      CALL Agrif_Update_Variable( trn_id, locupdate=(/1,2/), procname=updateTRC )
52# endif
53      !
54      Agrif_UseSpecialValueInUpdate = .FALSE.
55      !
56#endif
57      !
58   END SUBROUTINE Agrif_Update_Trc
59
60#ifdef key_vertical
61   SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
62      !!---------------------------------------------
63      !!           *** ROUTINE updateT ***
64      !!---------------------------------------------
65      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
66      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
67      LOGICAL, INTENT(in) :: before
68      !!
69      INTEGER :: ji,jj,jk,jn
70      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
71      REAL(wp) :: h_in(k1:k2)
72      REAL(wp) :: h_out(1:jpk)
73      INTEGER  :: N_in, N_out
74      REAL(wp) :: h_diff
75      REAL(wp) :: zrho_xy
76      REAL(wp) :: tabin(k1:k2,n1:n2)
77      !!---------------------------------------------
78      !
79      IF (before) THEN
80         AGRIF_SpecialValue = -999._wp
81         zrho_xy = Agrif_rhox() * Agrif_rhoy() 
82         DO jn = n1,n2-1
83            DO jk=k1,k2
84               DO jj=j1,j2
85                  DO ji=i1,i2
86                     tabres(ji,jj,jk,jn) = (    tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Kmm) ) &
87                       &                   * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
88                  END DO
89               END DO
90            END DO
91         END DO
92         DO jk=k1,k2
93            DO jj=j1,j2
94               DO ji=i1,i2
95                  tabres(ji,jj,jk,n2) =   tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) &
96                    &                  + (tmask(ji,jj,jk)-1)*999._wp
97               END DO
98            END DO
99         END DO
100      ELSE
101         tabres_child(:,:,:,:) = 0.
102         AGRIF_SpecialValue = 0._wp
103         DO jj=j1,j2
104            DO ji=i1,i2
105               N_in = 0
106               DO jk=k1,k2 !k2 = jpk of child grid
107                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT
108                  N_in = N_in + 1
109                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
110                  h_in(N_in) = tabres(ji,jj,jk,n2)
111               ENDDO
112               N_out = 0
113               DO jk=1,jpk ! jpk of parent grid
114                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
115                  N_out = N_out + 1
116                  h_out(N_out) = e3t(ji,jj,jk,Kmm) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
117               ENDDO
118               IF (N_in > 0) THEN !Remove this?
119                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
120                  IF (h_diff < -1.e-4) THEN
121                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
122                     print *,h_in(1:N_in)
123                     print *,h_out(1:N_out)
124                     STOP
125                  ENDIF
126                  DO jn=1,jptra
127                     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)
128                  ENDDO
129               ENDIF
130            ENDDO
131         ENDDO
132
133         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
134            ! Add asselin part
135            DO jn = 1,jptra
136               DO jk=1,jpk
137                  DO jj=j1,j2
138                     DO ji=i1,i2
139                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
140                           tr(ji,jj,jk,jn,Kbb) = ts(ji,jj,jk,jn,Kbb)                     & 
141                             &                  + atfp * ( tabres_child(ji,jj,jk,jn)     &
142                             &                                     - tr(ji,jj,jk,jn,Kmm) &
143                             &                           ) *      tmask(ji,jj,jk)
144                        ENDIF
145                     ENDDO
146                  ENDDO
147               ENDDO
148            ENDDO
149         ENDIF
150         DO jn = 1,jptra
151            DO jk=1,jpk
152               DO jj=j1,j2
153                  DO ji=i1,i2
154                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
155                        tr(ji,jj,jk,jn,Kmm) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
156                     END IF
157                  END DO
158               END DO
159            END DO
160         END DO
161      ENDIF
162      !
163   END SUBROUTINE updateTRC
164
165
166#else
167   SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
168      !!----------------------------------------------------------------------
169      !!                      *** ROUTINE updateTRC ***
170      !!----------------------------------------------------------------------
171      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
172      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
173      LOGICAL                                    , INTENT(in   ) ::   before
174      !!
175      INTEGER :: ji,jj,jk,jn
176      REAL(wp) :: ztb, ztnu, ztno
177      !!----------------------------------------------------------------------
178      !
179      !
180      IF (before) THEN
181         DO jn = n1,n2
182            DO jk=k1,k2
183               DO jj=j1,j2
184                  DO ji=i1,i2
185!> jc tmp
186                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Kmm) / e3t_0(ji,jj,jk)
187!                    tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Kmm)
188!< jc tmp
189                  END DO
190               END DO
191            END DO
192         END DO
193      ELSE
194!> jc tmp
195         DO jn = n1,n2
196            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
197              &                            * tmask(i1:i2,j1:j2,k1:k2)
198         ENDDO
199!< jc tmp
200         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
201            ! Add asselin part
202            DO jn = n1,n2
203               DO jk=k1,k2
204                  DO jj=j1,j2
205                     DO ji=i1,i2
206                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
207                           ztb  =     tr(ji,jj,jk,jn,Kbb) * e3t(ji,jj,jk,Kbb) ! fse3t_b prior update should be used
208                           ztnu = tabres(ji,jj,jk,jn)
209                           ztno =     tr(ji,jj,jk,jn,Kmm) * e3t(ji,jj,jk,Krhs)
210                           tr(ji,jj,jk,jn,Kbb) = ( ztb + atfp * ( ztnu - ztno) ) * tmask(ji,jj,jk)     &
211                             &                                                   /   e3t(ji,jj,jk,Kbb)
212                        ENDIF
213                     ENDDO
214                  ENDDO
215               ENDDO
216            ENDDO
217         ENDIF
218         DO jn = n1,n2
219            DO jk=k1,k2
220               DO jj=j1,j2
221                  DO ji=i1,i2
222                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
223                        tr(ji,jj,jk,jn,Kmm) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm)
224                     END IF
225                  END DO
226               END DO
227            END DO
228         END DO
229         !
230         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN
231            tr(i1:i2,j1:j2,k1:k2,n1:n2,Kbb)  = tr(i1:i2,j1:j2,k1:k2,n1:n2,Kmm)
232         ENDIF
233         !
234      ENDIF
235      !
236   END SUBROUTINE updateTRC
237#endif
238
239#else
240   !!----------------------------------------------------------------------
241   !!   Empty module                                           no TOP AGRIF
242   !!----------------------------------------------------------------------
243CONTAINS
244   SUBROUTINE agrif_top_update_empty
245      WRITE(*,*)  'agrif_top_update : You should not have seen this print! error?'
246   END SUBROUTINE agrif_top_update_empty
247#endif
248
249   !!======================================================================
250END MODULE agrif_top_update
Note: See TracBrowser for help on using the repository browser.