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_top_update.F90 in branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90 @ 9005

Last change on this file since 9005 was 9005, checked in by timgraham, 6 years ago

Add TOP update code

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