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.
vremap.F90 in NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90 @ 11741

Last change on this file since 11741 was 11741, checked in by jchanut, 5 years ago

#2222: correct definition of parent vertical grid on the child domain to perform vertical interpolation at boundaries. Use additionnal parent depths and number of levels arrays interpolated on the child grid domain to do so.
Correction of vertical interpolation of viscosity remains to be done as well as duplication of changes for passive tracers.

File size: 13.4 KB
Line 
1#undef PPR_LIB      /* USE PPR library */
2MODULE vremap
3!$AGRIF_DO_NOT_TREAT
4   !!======================================================================
5   !!                       ***  MODULE  vremap  ***
6   !! Ocean physics:  Vertical remapping routines
7   !!
8   !!======================================================================
9   !! History : 4.0  !  2019-09  (Jérôme Chanut)  Original code
10   !!----------------------------------------------------------------------
11   !!----------------------------------------------------------------------
12   !!
13   !!----------------------------------------------------------------------
14   USE par_oce
15#if defined PPR_LIB
16   USE ppr_1d   ! D. Engwirda piecewise polynomial reconstruction library
17#endif
18
19   IMPLICIT NONE
20   PRIVATE
21
22   PUBLIC   reconstructandremap
23
24   !!----------------------------------------------------------------------
25   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
26   !! $Id: vremap 11573 2019-09-19 09:18:03Z jchanut $
27   !! Software governed by the CeCILL license (see ./LICENSE)
28   !!----------------------------------------------------------------------
29CONTAINS
30
31#if ! defined PPR_LIB
32   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, kjpk_in, kjpk_out, kn_var)     
33      !!----------------------------------------------------------------------
34      !!                    ***  ROUTINE  reconstructandremap ***
35      !!
36      !! ** Purpose :   Brief description of the routine
37      !!
38      !! ** Method  :   description of the methodoloy used to achieve the
39      !!                objectives of the routine. Be as clear as possible!
40      !!
41      !! ** Action  : - first action (share memory array/varible modified
42      !!                in this routine
43      !!              - second action .....
44      !!              - .....
45      !!
46      !! References :   Author et al., Short_name_review, Year
47      !!                Give references if exist otherwise suppress these lines
48      !!-----------------------------------------------------------------------
49      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
50      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
51      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
52      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   phin       ! Input thicknesses
53      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   phout      ! Output thicknesses
54      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
55      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Remapped data
56      !
57      INTEGER             :: jk, jn, k1, kbox, ktop, ka, kbot
58      REAL(wp), PARAMETER :: dpthin = 1.D-3, dsmll = 1.0D-8
59      REAL(wp)            :: q, q01, q02, q001, q002, q0
60      REAL(wp)            :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
61      REAL(wp)            :: coeffremap(kjpk_in,3), zwork(kjpk_in,3), zwork2(kjpk_in+1,3)
62      REAL(wp)            :: z_win(1:kjpk_in+1), z_wout(1:kjpk_out+1)
63      !!-----------------------------------------------------------------------
64
65      z_win(1)=0._wp ; z_wout(1)= 0._wp
66      DO jk = 1, kjpk_in
67         z_win(jk+1)=z_win(jk)+phin(jk)
68      END DO
69     
70      DO jk = 1, kjpk_out
71         z_wout(jk+1)=z_wout(jk)+phout(jk)       
72      END DO       
73
74      DO jk = 2, kjpk_in
75         zwork(jk,1)=1._wp/(phin(jk-1)+phin(jk))
76      END DO
77       
78      DO jk = 2, kjpk_in-1
79         q0 = 1._wp / (phin(jk-1)+phin(jk)+phin(jk+1))
80         zwork(jk,2) = phin(jk-1) + 2._wp*phin(jk) + phin(jk+1)
81         zwork(jk,3) = q0
82      END DO       
83     
84      DO jn = 1, kn_var
85
86         DO jk = 2,kjpk_in
87            zwork2(jk,1) = zwork(jk,1)*(ptin(jk,jn)-ptin(jk-1,jn))
88         END DO
89       
90         coeffremap(:,1) = ptin(:,jn)
91 
92         DO jk = 2, kjpk_in-1
93            q001 = phin(jk)*zwork2(jk+1,1)
94            q002 = phin(jk)*zwork2(jk,1)       
95            IF (q001*q002 < 0._wp) then
96               q001 = 0._wp
97               q002 = 0._wp
98            ENDIF
99            q=zwork(jk,2)
100            q01=q*zwork2(jk+1,1)
101            q02=q*zwork2(jk,1)
102            IF (abs(q001) > abs(q02)) q001 = q02
103            IF (abs(q002) > abs(q01)) q002 = q01
104       
105            q=(q001-q002)*zwork(jk,3)
106            q001=q001-q*phin(jk+1)
107            q002=q002+q*phin(jk-1)
108       
109            coeffremap(jk,3)=coeffremap(jk,1)+q001
110            coeffremap(jk,2)=coeffremap(jk,1)-q002
111       
112            zwork2(jk,1)=(2._wp*q001-q002)**2
113            zwork2(jk,2)=(2._wp*q002-q001)**2
114         ENDDO
115       
116         DO jk = 1, kjpk_in
117            IF(jk.EQ.1 .OR. jk.EQ.kjpk_in .OR. phin(jk).LE.dpthin) THEN
118               coeffremap(jk,3) = coeffremap(jk,1)
119               coeffremap(jk,2) = coeffremap(jk,1)
120               zwork2(jk,1) = 0._wp
121               zwork2(jk,2) = 0._wp
122            ENDIF
123         END DO
124       
125         DO jk = 2, kjpk_in
126            q002 = max(zwork2(jk-1,2),dsmll)
127            q001 = max(zwork2(jk,1)  ,dsmll)
128            zwork2(jk,3) = (q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)
129         END DO
130       
131         zwork2(1,3) = 2._wp*coeffremap(1,1)-zwork2(2,3)
132         zwork2(kjpk_in+1,3)=2._wp*coeffremap(kjpk_in,1)-zwork2(kjpk_in,3)
133 
134         DO jk = 1, kjpk_in
135            q01=zwork2(jk+1,3)-coeffremap(jk,1)
136            q02=coeffremap(jk,1)-zwork2(jk,3)
137            q001=2._wp*q01
138            q002=2._wp*q02
139            IF (q01*q02<0._wp) then
140               q01=0._wp
141               q02=0._wp
142            ELSEIF (abs(q01)>abs(q002)) then
143               q01=q002
144            ELSEIF (abs(q02)>abs(q001)) then
145               q02=q001
146            ENDIF
147            coeffremap(jk,2)=coeffremap(jk,1)-q02
148            coeffremap(jk,3)=coeffremap(jk,1)+q01
149         ENDDO
150
151         zbot=0._wp
152         kbot=1
153         DO jk=1,kjpk_out
154            ztop=zbot  !top is bottom of previous layer
155            ktop=kbot
156            IF     (ztop.GE.z_win(ktop+1)) then
157               ktop=ktop+1
158            ENDIF
159       
160            zbot=z_wout(jk+1)
161            zthk=zbot-ztop
162
163            IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(kjpk_out+1)) THEN
164 
165               kbot=ktop
166               DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.kjpk_in)
167                  kbot=kbot+1
168               ENDDO
169               zbox=zbot
170               DO k1= jk+1,kjpk_out
171                  IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN
172                     exit !thick layer
173                  ELSE
174                     zbox=z_wout(k1+1)  !include thin adjacent layers
175                     IF(zbox.EQ.z_wout(kjpk_out+1)) THEN
176                        exit !at bottom
177                     ENDIF
178                  ENDIF
179               ENDDO
180               zthk=zbox-ztop
181
182               kbox=ktop
183               DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.kjpk_in)
184                  kbox=kbox+1
185               ENDDO
186         
187               IF(ktop.EQ.kbox) THEN
188                  IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN
189                     IF(phin(kbox).GT.dpthin) THEN
190                        q001 = (zbox-z_win(kbox))/phin(kbox)
191                        q002 = (ztop-z_win(kbox))/phin(kbox)
192                        q01=q001**2+q002**2+q001*q002+1._wp-2._wp*(q001+q002)
193                        q02=q01-1._wp+(q001+q002)
194                        q0=1._wp-q01-q02
195                     ELSE
196                        q0  = 1._wp
197                        q01 = 0._wp
198                        q02 = 0._wp
199                     ENDIF
200                     ptout(jk,jn)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
201                  ELSE
202                     ptout(jk,jn) = ptin(kbox,jn)
203                  ENDIF
204               ELSE
205                  IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
206                     ka = jk
207                  ELSEIF (kbox-ktop.GE.3) THEN
208                     ka = (kbox+ktop)/2
209                  ELSEIF (phin(ktop).GE.phin(kbox)) THEN
210                     ka = ktop
211                  ELSE
212                     ka = kbox
213                  ENDIF !choose ka
214   
215                  offset=coeffremap(ka,1)
216   
217                  qtop = z_win(ktop+1)-ztop !partial layer thickness
218                  IF(phin(ktop).GT.dpthin) THEN
219                     q=(ztop-z_win(ktop))/phin(ktop)
220                     q01=q*(q-1._wp)
221                     q02=q01+q
222                     q0=1._wp-q01-q02           
223                  ELSE
224                     q0  = 1._wp
225                     q01 = 0._wp
226                     q02 = 0._wp
227                  ENDIF
228               
229                  tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
230   
231                  DO k1= ktop+1,kbox-1
232                     tsum =tsum +(coeffremap(k1,1)-offset)*phin(k1)
233                  ENDDO !k1
234               
235                  qbot = zbox-z_win(kbox) !partial layer thickness
236                  IF(phin(kbox).GT.dpthin) THEN
237                     q=qbot/phin(kbox)
238                     q01=(q-1._wp)**2
239                     q02=q01-1._wp+q
240                     q0=1_wp-q01-q02                           
241                  ELSE
242                     q0  = 1._wp
243                     q01 = 0._wp
244                     q02 = 0._wp
245                  ENDIF
246             
247                  tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
248               
249                  rpsum=1._wp / zthk
250                  ptout(jk,jn)=offset+tsum*rpsum
251                 
252               ENDIF !single or multiple layers
253            ELSE
254               IF (jk==1) THEN
255                  write(*,'(a7,i4,i4,3f12.5)')'problem = ',kjpk_in,kjpk_out,zthk,z_wout(jk+1),phout(1)
256               ENDIF
257               ptout(jk,jn) = ptout(jk-1,jn)
258             
259            ENDIF !normal:thin layer
260         ENDDO !jk
261
262      END DO ! loop over variables
263           
264   END SUBROUTINE reconstructandremap
265
266#else
267
268   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, kjpk_in, kjpk_out, kn_var)
269      !!----------------------------------------------------------------------
270      !!                    *** ROUTINE  reconstructandremap ***
271      !!
272      !! ** Purpose :   Conservative remapping of a vertical column
273      !!                from one set of layers to an other one.
274      !!
275      !! ** Method  :   Uses D. Engwirda Piecewise Polynomial Reconstruction library.
276      !!                https://github.com/dengwirda/PPR
277      !!               
278      !!
279      !! References :   Engwirda, Darren & Kelley, Maxwell. (2015). A WENO-type
280      !!                slope-limiter for a family of piecewise polynomial methods.
281      !!                https://arxiv.org/abs/1606.08188
282      !!-----------------------------------------------------------------------
283      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
284      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
285      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
286      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   phin       ! Input thicknesses
287      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   phout      ! Output thicknesses
288      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
289      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Remapped data
290      !
291      INTEGER, PARAMETER :: ndof = 1
292      INTEGER  :: jk, jn
293      REAL(wp) ::  zwin(kjpk_in+1) ,  ztin(ndof, kn_var, kjpk_in)
294      REAL(wp) :: zwout(kjpk_out+1), ztout(ndof, kn_var, kjpk_out)
295      TYPE(rmap_work) :: work
296      TYPE(rmap_opts) :: opts
297      TYPE(rcon_ends) :: bc_l(kn_var)
298      TYPE(rcon_ends) :: bc_r(kn_var)
299      !!--------------------------------------------------------------------
300     
301      ! Set interfaces and input data:
302      zwin(1) = 0._wp
303      DO jk = 2, kjpk_in + 1
304         zwin(jk) = zwin(jk-1) + phin(jk-1) 
305      END DO
306     
307      DO jn = 1, kn_var 
308         DO jk = 1, kjpk_in
309            ztin(ndof, jn, jk) =  ptin(jk, jn)
310         END DO
311      END DO
312
313      zwout(1) = 0._wp
314      DO jk = 2, kjpk_out + 1
315         zwout(jk) = zwout(jk-1) + phout(jk-1) 
316      END DO
317
318      ! specify methods
319!      opts%edge_meth = p1e_method     ! 1st-order edge interp.
320!      opts%cell_meth = plm_method     ! PLM method in cells
321      opts%edge_meth = p3e_method     ! 3rd-order edge interp.
322      opts%cell_meth = ppm_method     ! PPM method in cells   
323!      opts%edge_meth = p5e_method     ! 5th-order edge interp.
324!      opts%cell_meth = pqm_method     ! PQM method in cells
325
326      ! limiter
327!      opts%cell_lims = null_limit     ! no lim.
328      opts%cell_lims = mono_limit     ! monotone limiter   
329 
330      ! set boundary conditions
331      bc_l%bcopt = bcon_loose         ! "loose" = extrapolate
332      bc_r%bcopt = bcon_loose
333!      bc_l%bcopt = bcon_slope       
334!      bc_r%bcopt = bcon_slope
335
336      ! init. method workspace
337      CALL work%init(kjpk_in+1, kn_var, opts)
338
339      ! remap
340      CALL rmap1d(kjpk_in+1, kjpk_out+1, kn_var, ndof, &
341      &           zwin, zwout, ztin, ztout,            &
342      &           bc_l, bc_r, work, opts)
343
344      ! clear method workspace
345      CALL work%free()
346
347      DO jn = 1, kn_var 
348         DO jk = 1, kjpk_out
349            ptout(jk, jn) = ztout(1, jn, jk)
350         END DO
351      END DO
352           
353   END SUBROUTINE reconstructandremap
354#endif
355
356   !!======================================================================
357!$AGRIF_END_DO_NOT_TREAT
358END MODULE vremap
Note: See TracBrowser for help on using the repository browser.