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

Last change on this file since 13226 was 12377, checked in by acc, 12 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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   !! * 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 = plm_method     ! PLM method in cells
323      opts%edge_meth = p3e_method     ! 3rd-order edge interp.
324      opts%cell_meth = ppm_method     ! PPM method in cells   
325!      opts%edge_meth = p5e_method     ! 5th-order edge interp.
326!      opts%cell_meth = pqm_method     ! PQM method in cells
327
328      ! limiter
329!      opts%cell_lims = null_limit     ! no lim.
330      opts%cell_lims = mono_limit     ! monotone limiter   
331 
332      ! set boundary conditions
333      bc_l%bcopt = bcon_loose         ! "loose" = extrapolate
334      bc_r%bcopt = bcon_loose
335!      bc_l%bcopt = bcon_slope       
336!      bc_r%bcopt = bcon_slope
337
338      ! init. method workspace
339      CALL work%init(kjpk_in+1, kn_var, opts)
340
341      ! remap
342      CALL rmap1d(kjpk_in+1, kjpk_out+1, kn_var, ndof, &
343      &           zwin, zwout, ztin, ztout,            &
344      &           bc_l, bc_r, work, opts)
345
346      ! clear method workspace
347      CALL work%free()
348
349      DO jn = 1, kn_var 
350         DO jk = 1, kjpk_out
351            ptout(jk, jn) = ztout(1, jn, jk)
352         END DO
353      END DO
354           
355   END SUBROUTINE reconstructandremap
356#endif
357
358   SUBROUTINE remap_linear(ptin, pzin, ptout, pzout, kjpk_in, kjpk_out, kn_var)
359      !!----------------------------------------------------------------------
360      !!                    *** ROUTINE  remap_linear ***
361      !!
362      !! ** Purpose :   Linear interpolation based on input/ouputs depths
363      !!
364      !!-----------------------------------------------------------------------
365      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
366      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
367      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
368      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   pzin       ! Input depths
369      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   pzout      ! Output depths
370      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
371      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Interpolated data
372      !
373      INTEGER  :: jkin, jkout, jn
374      !!--------------------------------------------------------------------
375      !     
376      DO jkout = 1, kjpk_out !  Loop over destination grid
377         !
378         IF     ( pzout(jkout) <=  pzin(  1    ) ) THEN  ; ptout(jkout,1:kn_var) = ptin(    1    ,1:kn_var)         
379         ELSEIF ( pzout(jkout) >= pzin(kjpk_in) ) THEN   ; ptout(jkout,1:kn_var) = ptin( kjpk_in ,1:kn_var)
380         ELSEIF ( ( pzout(jkout) > pzin(1) ).AND.( pzout(jkout) < pzin(kjpk_in) )) THEN
381            DO jkin = 1, kjpk_in - 1 !  Loop over source grid
382               IF ( pzout(jkout) < pzin(jkin+1) ) THEN
383                  DO jn = 1, kn_var
384                     ptout(jkout,jn) =  ptin(jkin,jn) + &
385                                     & (pzout(jkout) - pzin(jkin)) / (pzin(jkin+1)    - pzin(jkin)) &
386                                     &                             * (ptin(jkin+1,jn) - ptin(jkin,jn))
387                  END DO 
388                  EXIT
389               ENDIF 
390            END DO
391         ENDIF
392         !
393      END DO
394
395   END SUBROUTINE remap_linear
396
397   !!======================================================================
398!$AGRIF_END_DO_NOT_TREAT
399END MODULE vremap
Note: See TracBrowser for help on using the repository browser.