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

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

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

#2222, turn on the use of PPR library with vertical nesting

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