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

source: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/NST/vremap.F90 @ 13463

Last change on this file since 13463 was 13463, checked in by andmirek, 4 years ago

Ticket #2195:update to trunk 13461
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.