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_interp.F90 in NEMO/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_top_interp.F90

Last change on this file was 14218, checked in by jchanut, 3 years ago

#2222, Fixes uninitialized arrays with vertical remap
  • Property svn:keywords set to Id
File size: 8.7 KB
Line 
1MODULE agrif_top_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_top_interp  ***
4   !! AGRIF: interpolation package for TOP
5   !!======================================================================
6   !! History :  2.0  !  ???
7   !!----------------------------------------------------------------------
8#if defined key_agrif && defined key_top
9   !!----------------------------------------------------------------------
10   !!   'key_agrif'                                              AGRIF zoom
11   !!   'key_top'                                           on-line tracers
12   !!----------------------------------------------------------------------
13   USE par_oce
14   USE oce
15   USE dom_oce     
16   USE agrif_oce
17   USE agrif_top_sponge
18   USE par_trc
19   USE trc
20   USE vremap
21   !
22   USE lib_mpp     ! MPP library
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC Agrif_trc, interptrn
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31  !!----------------------------------------------------------------------
32   !! NEMO/NST 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE Agrif_trc
39      !!----------------------------------------------------------------------
40      !!                   ***  ROUTINE Agrif_trc  ***
41      !!----------------------------------------------------------------------
42      !
43      IF( Agrif_Root() )   RETURN
44      !
45      Agrif_SpecialValue    = 0._wp
46      Agrif_UseSpecialValue = .TRUE.
47      l_vremap              = ln_vert_remap
48      !
49      CALL Agrif_Bc_variable( trn_id, procname=interptrn )
50      !
51      Agrif_UseSpecialValue = .FALSE.
52      l_vremap              = .FALSE.
53      !
54   END SUBROUTINE Agrif_trc
55
56   SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
57      !!----------------------------------------------------------------------
58      !!                  *** ROUTINE interptrn ***
59      !!----------------------------------------------------------------------
60      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
61      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
62      LOGICAL                                     , INTENT(in   ) ::   before
63      !
64      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
65      INTEGER  ::   N_in, N_out
66      INTEGER  :: item
67      ! vertical interpolation:
68      REAL(wp) :: zhtot, zwgt
69      REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin, tabin_i
70      REAL(wp), DIMENSION(k1:k2) :: z_in, h_in_i, z_in_i
71      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
72      !!----------------------------------------------------------------------
73
74      IF( before ) THEN
75
76         item = Kmm_a
77         IF( l_ini_child )   Kmm_a = Kbb_a 
78
79         DO jn = 1,jptra
80            DO jk=k1,k2
81               DO jj=j1,j2
82                 DO ji=i1,i2
83                       ptab(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)
84                 END DO
85              END DO
86           END DO
87         END DO
88
89         IF( l_vremap .OR. l_ini_child .OR. ln_zps ) THEN
90            ! Fill cell depths (i.e. gdept) to be interpolated
91            ! Warning: these are masked, hence extrapolated prior interpolation.
92            DO jj=j1,j2
93               DO ji=i1,i2
94                  ptab(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kmm_a)
95                  DO jk=k1+1,k2
96                     ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * &
97                        & ( ptab(ji,jj,jk-1,jptra+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kmm_a)+e3t(ji,jj,jk,Kmm_a)) )
98                  END DO
99               END DO
100            END DO
101
102            ! Save ssh at last level:
103            IF (.NOT.ln_linssh) THEN
104               ptab(i1:i2,j1:j2,k2,jptra+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1) 
105            END IF     
106         ENDIF
107         Kmm_a = item
108
109      ELSE
110         item = Krhs_a
111         IF( l_ini_child )   Krhs_a = Kbb_a 
112
113         IF( l_vremap .OR. l_ini_child ) THEN
114            IF (ln_linssh) ptab(i1:i2,j1:j2,k2,n2) = 0._wp 
115               
116            DO jj=j1,j2
117               DO ji=i1,i2
118                  tr(ji,jj,:,:,Krhs_a) = 0. 
119                  !
120                  ! Build vertical grids:
121                  N_in = mbkt_parent(ji,jj)
122                  N_out = mbkt(ji,jj)
123                  IF (N_in*N_out > 0) THEN
124                     ! Input grid (account for partial cells if any):
125                     DO jk=1,N_in
126                        z_in(jk) = ptab(ji,jj,jk,n2) - ptab(ji,jj,k2,n2)
127                        tabin(jk,1:jptra) = ptab(ji,jj,jk,1:jptra)
128                     END DO
129                 
130                     ! Intermediate grid:
131                     IF ( l_vremap ) THEN
132                        DO jk = 1, N_in
133                           h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
134                             &       (1._wp + ptab(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
135                        END DO
136                        z_in_i(1) = 0.5_wp * h_in_i(1)
137                        DO jk=2,N_in
138                           z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
139                        END DO
140                        z_in_i(1:N_in) = z_in_i(1:N_in)  - ptab(ji,jj,k2,n2)
141                     ENDIF
142
143                     ! Output (Child) grid:
144                     DO jk=1,N_out
145                        h_out(jk) = e3t(ji,jj,jk,Krhs_a)
146                     END DO
147                     z_out(1) = 0.5_wp * h_out(1)
148                     DO jk=2,N_out
149                        z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
150                     END DO
151                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a)               
152
153                     IF( l_ini_child ) THEN
154                        CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a),          &
155                                      &   z_out(1:N_out),N_in,N_out,jptra) 
156                     ELSE 
157                        CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabin_i(1:N_in,1:jptra),                     &
158                                     &   z_in_i(1:N_in),N_in,N_in,jptra)
159                        CALL reconstructandremap(tabin_i(1:N_in,1:jptra),h_in_i(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a), &
160                                      &   h_out(1:N_out),N_in,N_out,jptra)   
161                     ENDIF
162                  ENDIF
163               END DO
164            END DO
165            Krhs_a = item
166 
167         ELSE
168
169            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
170                                             ! linear vertical interpolation
171               DO jj=j1,j2
172                  DO ji=i1,i2
173                     !
174                     N_in  = mbkt(ji,jj)
175                     N_out = mbkt(ji,jj)
176                     z_in(1) = ptab(ji,jj,1,n2)
177                     tabin(1,1:jptra) = ptab(ji,jj,1,1:jptra)
178                     DO jk=2, N_in
179                        z_in(jk) = ptab(ji,jj,jk,n2)
180                        tabin(jk,1:jptra) = ptab(ji,jj,jk,1:jptra)
181                     END DO
182                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2)
183                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a)
184                     DO jk=2, N_out
185                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a))
186                     END DO
187                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a)
188                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jptra), &
189                                   &   z_out(1:N_out),N_in,N_out,jptra) 
190                  END DO
191               END DO
192
193            ENDIF
194
195            DO jn=1, jptra
196                tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
197            END DO
198         ENDIF
199
200      ENDIF
201      !
202   END SUBROUTINE interptrn
203
204#else
205   !!----------------------------------------------------------------------
206   !!   Empty module                                           no TOP AGRIF
207   !!----------------------------------------------------------------------
208CONTAINS
209   SUBROUTINE Agrif_TOP_Interp_empty
210      !!---------------------------------------------
211      !!   *** ROUTINE agrif_Top_Interp_empty ***
212      !!---------------------------------------------
213      WRITE(*,*)  'agrif_top_interp : You should not have seen this print! error?'
214   END SUBROUTINE Agrif_TOP_Interp_empty
215#endif
216
217   !!======================================================================
218END MODULE agrif_top_interp
Note: See TracBrowser for help on using the repository browser.