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

Last change on this file since 11802 was 11802, checked in by jchanut, 9 months ago

#2222, 1) add linear interpolation in vremap module.
2) Switch remapping of viscosity from polynomial to linear.
3) Move to truly volume weighted averages for parent to child update.

File size: 15.5 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, remap_linear
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   SUBROUTINE remap_linear(ptin, pzin, ptout, pzout, kjpk_in, kjpk_out, kn_var)
357      !!----------------------------------------------------------------------
358      !!                    *** ROUTINE  remap_linear ***
359      !!
360      !! ** Purpose :   Linear interpolation based on input/ouputs depths
361      !!
362      !!-----------------------------------------------------------------------
363      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
364      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
365      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
366      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   pzin       ! Input depths
367      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   pzout      ! Output depths
368      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
369      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Interpolated data
370      !
371      INTEGER  :: jkin, jkout, jn
372      !!--------------------------------------------------------------------
373      !     
374      DO jkout = 1, kjpk_out !  Loop over destination grid
375         !
376         IF     ( pzout(jkout) < pzin(  1    ) ) THEN    ; ptout(jkout,1:kn_var) = ptin(    1    ,1:kn_var)         
377         ELSEIF ( pzout(jkout) > pzin(kjpk_in) ) THEN    ; ptout(jkout,1:kn_var) = ptin( kjpk_in ,1:kn_var)
378         ELSEIF ( ( pzout(jkout) >= pzin(1) ).AND.( pzout(jkout) <= pzin(kjpk_in) )) THEN
379            DO jkin = 1, kjpk_in - 1 !  Loop over source grid
380               IF ( pzout(jkout) >=  pzin(jkin) ) THEN
381                  DO jn = 1, kn_var
382                     ptout(jkout,jn) =  ptin(jkin,jn) + &
383                                     & (pzout(jkout) - pzin(jkin)) / (pzin(jkin+1)    - pzin(jkin)) &
384                                     &                             * (ptin(jkin+1,jn) - ptin(jkin,jn))
385                  END DO 
386                  EXIT
387               ENDIF 
388            END DO
389         ENDIF
390         !
391      END DO
392           
393   END SUBROUTINE remap_linear
394
395   !!======================================================================
396!$AGRIF_END_DO_NOT_TREAT
397END MODULE vremap
Note: See TracBrowser for help on using the repository browser.