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 NEMO/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_top_update.F90 @ 12808

Last change on this file since 12808 was 12489, checked in by davestorkey, 4 years ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp -> rn_atfp (namelist parameter used everywhere) and rau0 -> rho0.

This version gives bit-comparable results to the previous version of the trunk.
  • Property svn:keywords set to Id
File size: 9.5 KB
RevLine 
[5656]1#undef DECAL_FEEDBACK
[628]2
[636]3MODULE agrif_top_update
[6140]4   !!======================================================================
5   !!                ***  MODULE agrif_top_update  ***
[9019]6   !! AGRIF :   update package for passive tracers (TOP)
7   !!======================================================================
[6140]8   !! History : 
9   !!----------------------------------------------------------------------
[1206]10#if defined key_agrif && defined key_top
[9019]11   !!----------------------------------------------------------------------
12   !!   'key_agrif'                                              AGRIF zoom
13   !!   'key_TOP'                                           on-line tracers
14   !!----------------------------------------------------------------------
[636]15   USE par_oce
16   USE oce
[9019]17   USE dom_oce
18   USE agrif_oce
[6140]19   USE par_trc
20   USE trc
[12377]21   USE vremap
[628]22
[636]23   IMPLICIT NONE
24   PRIVATE
[628]25
[636]26   PUBLIC Agrif_Update_Trc
[628]27
[1156]28   !!----------------------------------------------------------------------
[9598]29   !! NEMO/NST 4.0 , NEMO Consortium (2018)
[1156]30   !! $Id$
[10068]31   !! Software governed by the CeCILL license (see ./LICENSE)
[1156]32   !!----------------------------------------------------------------------
[5656]33CONTAINS
[636]34
[9031]35   SUBROUTINE Agrif_Update_Trc( )
[6140]36      !!----------------------------------------------------------------------
37      !!                   *** ROUTINE Agrif_Update_Trc ***
38      !!----------------------------------------------------------------------
[5656]39      !
[9031]40      IF (Agrif_Root()) RETURN 
41      !
[628]42      Agrif_UseSpecialValueInUpdate = .TRUE.
[6140]43      Agrif_SpecialValueFineGrid    = 0._wp
[5656]44      !
45# if ! defined DECAL_FEEDBACK
[9031]46      CALL Agrif_Update_Variable(trn_id, procname=updateTRC )
47!      CALL Agrif_Update_Variable( trn_id, locupdate=(/0,2/), procname=updateTRC )
[5656]48# else
[9031]49      CALL Agrif_Update_Variable(trn_id, locupdate=(/1,0/),procname=updateTRC )
50!      CALL Agrif_Update_Variable( trn_id, locupdate=(/1,2/), procname=updateTRC )
[5656]51# endif
52      !
[628]53      Agrif_UseSpecialValueInUpdate = .FALSE.
[9031]54      !
[636]55   END SUBROUTINE Agrif_Update_Trc
[628]56
[9031]57#ifdef key_vertical
58   SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
59      !!---------------------------------------------
60      !!           *** ROUTINE updateT ***
61      !!---------------------------------------------
62      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
63      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
64      LOGICAL, INTENT(in) :: before
65      !!
66      INTEGER :: ji,jj,jk,jn
[12377]67      REAL(wp) :: ztb, ztnu, ztno
[9031]68      REAL(wp) :: h_in(k1:k2)
69      REAL(wp) :: h_out(1:jpk)
70      INTEGER  :: N_in, N_out
71      REAL(wp) :: h_diff
[12377]72      REAL(wp) :: tabin(k1:k2,1:jptra)
73      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jptra) :: tabres_child
[9031]74      !!---------------------------------------------
75      !
76      IF (before) THEN
77         AGRIF_SpecialValue = -999._wp
78         DO jn = n1,n2-1
79            DO jk=k1,k2
80               DO jj=j1,j2
81                  DO ji=i1,i2
[12377]82                     tabres(ji,jj,jk,jn) = (tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) &
[9031]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
[12377]91                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
[9031]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
[12377]112                  h_out(N_out) = e3t(ji,jj,jk,Kmm_a) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
[9031]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
[12377]122                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra)
[9031]123               ENDIF
124            ENDDO
125         ENDDO
[12377]126         !
[12489]127         IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
[9031]128            ! Add asselin part
129            DO jn = 1,jptra
[12377]130               DO jk=1,jpkm1
[9031]131                  DO jj=j1,j2
132                     DO ji=i1,i2
133                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
[12377]134                           ztb  = tr(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
135                           ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
136                           ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
[12489]137                           tr(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
[12377]138                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
[9031]139                        ENDIF
140                     ENDDO
141                  ENDDO
142               ENDDO
143            ENDDO
144         ENDIF
145         DO jn = 1,jptra
[12377]146            DO jk=1,jpkm1
[9031]147               DO jj=j1,j2
148                  DO ji=i1,i2
149                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
[12377]150                        tr(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn)
[9031]151                     END IF
152                  END DO
153               END DO
154            END DO
155         END DO
[12377]156         !
[12489]157         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]158            tr(i1:i2,j1:j2,1:jpkm1,1:jptra,Kbb_a)  = tr(i1:i2,j1:j2,1:jpkm1,1:jptra,Kmm_a)
159         ENDIF
160         !
161
[9031]162      ENDIF
163      !
164   END SUBROUTINE updateTRC
165
166
167#else
168   SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
[6140]169      !!----------------------------------------------------------------------
[9031]170      !!                      *** ROUTINE updateTRC ***
[6140]171      !!----------------------------------------------------------------------
172      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
[9031]173      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
[6140]174      LOGICAL                                    , INTENT(in   ) ::   before
[5656]175      !!
[9031]176      INTEGER :: ji,jj,jk,jn
177      REAL(wp) :: ztb, ztnu, ztno
[6140]178      !!----------------------------------------------------------------------
[5656]179      !
[9031]180      !
181      IF (before) THEN
182         DO jn = n1,n2
183            DO jk=k1,k2
184               DO jj=j1,j2
185                  DO ji=i1,i2
186!> jc tmp
[12377]187                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
188!                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a)
[9031]189!< jc tmp
190                  END DO
191               END DO
192            END DO
193         END DO
[5656]194      ELSE
[9031]195!> jc tmp
196         DO jn = n1,n2
197            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
198                                         & * tmask(i1:i2,j1:j2,k1:k2)
199         ENDDO
200!< jc tmp
[12489]201         IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
[4491]202            ! Add asselin part
[5656]203            DO jn = n1,n2
[9031]204               DO jk=k1,k2
205                  DO jj=j1,j2
206                     DO ji=i1,i2
207                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
[12377]208                           ztb  = tr(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
[9031]209                           ztnu = tabres(ji,jj,jk,jn)
[12377]210                           ztno = tr(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
[12489]211                           tr(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
[12377]212                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
[1300]213                        ENDIF
[9031]214                     ENDDO
215                  ENDDO
216               ENDDO
217            ENDDO
[636]218         ENDIF
[9031]219         DO jn = n1,n2
220            DO jk=k1,k2
221               DO jj=j1,j2
222                  DO ji=i1,i2
223                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
[12377]224                        tr(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
[5656]225                     END IF
226                  END DO
227               END DO
228            END DO
229         END DO
[9031]230         !
[12489]231         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
[12377]232            tr(i1:i2,j1:j2,k1:k2,n1:n2,Kbb_a)  = tr(i1:i2,j1:j2,k1:k2,n1:n2,Kmm_a)
[9031]233         ENDIF
234         !
[5656]235      ENDIF
236      !
[636]237   END SUBROUTINE updateTRC
[9031]238#endif
[628]239
240#else
[9019]241   !!----------------------------------------------------------------------
242   !!   Empty module                                           no TOP AGRIF
243   !!----------------------------------------------------------------------
[636]244CONTAINS
245   SUBROUTINE agrif_top_update_empty
246      WRITE(*,*)  'agrif_top_update : You should not have seen this print! error?'
247   END SUBROUTINE agrif_top_update_empty
[628]248#endif
[6140]249
250   !!======================================================================
[5656]251END MODULE agrif_top_update
Note: See TracBrowser for help on using the repository browser.