source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_update.F90 @ 11574

Last change on this file since 11574 was 11574, checked in by jchanut, 12 months ago

#2222, import changes from dev_r10973_AGRIF-01_jchanut_small_jpi_jpj (i.e. #2199)

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