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.
Changeset 11603 – NEMO

Changeset 11603


Ignore:
Timestamp:
2019-09-26T17:27:43+02:00 (5 years ago)
Author:
jchanut
Message:

#2222, 1) correct time interpolation of barotropic velocities in corners. 2) Clean remapping module and enable remapping several variables at the same time. At this stage, vertical remapping doesn't change VORTEX results with an identical vertical grid ONLY in one way mode and a linearized free surface (within truncature errors).

Location:
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90

    r11590 r11603  
    4545   PUBLIC   interpe3t, interpumsk, interpvmsk 
    4646 
     47   INTEGER ::   bdy_tinterp = 0 
     48 
    4749#  include "vectopt_loop_substitute.h90" 
    4850   !!---------------------------------------------------------------------- 
    … …  
    501503         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )  
    502504         ! 
     505         bdy_tinterp = 1 
    503506         CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After 
    504          CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  ) 
     507         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )   
    505508         ! 
     509         bdy_tinterp = 2 
    506510         CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before 
    507          CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )          
     511         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )    
    508512      ELSE ! Linear interpolation 
    509513         ! 
    … …  
    703707               ENDDO 
    704708               IF (N_in > 0) THEN 
    705                   DO jn=1,jpts 
    706                      call reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),ptab_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    707                   ENDDO 
     709                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ptab_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
    708710               ENDIF 
    709711            ENDDO 
    … …  
    800802              ENDIF 
    801803          
    802               IF (N_in * N_out > 0) THEN 
    803                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    804 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
    805                  if (h_diff < -1.e4) then 
    806                     print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
    807 !                    stop 
    808                  endif 
    809               ENDIF 
    810               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     804              CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    811805            ENDDO 
    812806         ENDDO 
    … …  
    881875                 CYCLE 
    882876               ENDIF 
    883                call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     877               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    884878            END DO 
    885879         END DO 
    … …  
    916910         DO ji = i1, i2 
    917911            DO jj = j1, j2 
    918                IF    ( utint_stage(ji,jj) == 1  ) THEN 
    919                   ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        & 
    920                      &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
    921                ELSEIF( utint_stage(ji,jj) == 2  ) THEN 
    922                   ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
    923                      &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
    924                ELSEIF( utint_stage(ji,jj) == 0  ) THEN                 
    925                   ztcoeff = 1._wp 
    926                ELSE 
    927                   ztcoeff = 0._wp 
    928                ENDIF 
    929                !    
    930                ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj) 
    931                !             
    932                IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN 
    933                   ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1) 
    934                   utint_stage(ji,jj) = 3 
    935                ELSE 
     912               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN 
     913                  IF    ( utint_stage(ji,jj) == 1  ) THEN 
     914                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        & 
     915                        &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
     916                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN 
     917                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
     918                        &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
     919                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN                 
     920                     ztcoeff = 1._wp 
     921                  ELSE 
     922                     ztcoeff = 0._wp 
     923                  ENDIF 
     924                  !    
     925                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj) 
     926                  !             
     927                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN 
     928                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1) 
     929                  ENDIF 
     930                  ! 
    936931                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1 
    937932               ENDIF 
    … …  
    966961         DO ji = i1, i2 
    967962            DO jj = j1, j2 
    968                IF    ( vtint_stage(ji,jj) == 1  ) THEN 
    969                   ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        & 
    970                      &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
    971                ELSEIF( vtint_stage(ji,jj) == 2  ) THEN 
    972                   ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
    973                      &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
    974                ELSEIF( vtint_stage(ji,jj) == 0  ) THEN                 
    975                   ztcoeff = 1._wp 
    976                ELSE 
    977                   ztcoeff = 0._wp 
    978                ENDIF 
    979                !    
    980                vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj) 
    981                !             
    982                IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN 
    983                   vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1) 
    984                   vtint_stage(ji,jj) = 3 
    985                ELSE 
     963               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN 
     964                  IF    ( vtint_stage(ji,jj) == 1  ) THEN 
     965                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        & 
     966                        &               - zt0**2._wp * (       zt0 - 1._wp)        ) 
     967                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN 
     968                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp & 
     969                        &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
     970                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN                 
     971                     ztcoeff = 1._wp 
     972                  ELSE 
     973                     ztcoeff = 0._wp 
     974                  ENDIF 
     975                  !    
     976                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj) 
     977                  !             
     978                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN 
     979                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1) 
     980                  ENDIF 
     981                  ! 
    986982                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1 
    987983               ENDIF 
    … …  
    12491245               ENDDO 
    12501246               IF (N_in > 0) THEN 
    1251                   CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out) 
     1247                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out,1) 
    12521248               ENDIF 
    12531249            ENDDO 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90

    r11590 r11603  
    346346               ENDDO 
    347347               IF (N_in > 0) THEN 
    348                   DO jn=1,jpts 
    349                      call reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    350                   ENDDO 
     348                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
    351349               ENDIF 
    352350            ENDDO 
    … …  
    488486                 endif 
    489487              ENDIF 
    490               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     488              CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    491489          
    492490            ENDDO 
    … …  
    635633                 endif 
    636634              ENDIF 
    637               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     635              CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    638636            ENDDO 
    639637         ENDDO 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90

    r11590 r11603  
    339339                     STOP 
    340340                  ENDIF 
    341                   DO jn=n1,n2-1 
    342                      CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    343                   ENDDO 
     341                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
    344342               ENDIF 
    345343            ENDDO 
    … …  
    524522                     ENDDO 
    525523                  ENDIF 
    526                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     524                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    527525                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
    528526               ENDIF 
    … …  
    709707                     ENDDO 
    710708                  ENDIF 
    711                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     709                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    712710                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
    713711               ENDIF 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_interp.F90

    r11590 r11603  
    106106               ENDDO 
    107107               IF (N_in > 0) THEN 
    108                   DO jn=1,jptra 
    109                      call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
    110                   ENDDO 
     108                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,ptab_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra) 
    111109               ENDIF 
    112110            ENDDO 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_sponge.F90

    r11590 r11603  
    117117               ENDDO 
    118118               IF (N_in > 0) THEN 
    119                   DO jn=1,jptra 
    120                      call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
    121                   ENDDO 
     119                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,tabres_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra) 
    122120               ENDIF 
    123121            ENDDO 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_update.F90

    r11574 r11603  
    120120                     STOP 
    121121                  ENDIF 
    122                   DO jn=1,jptra 
    123                      CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    124                   ENDDO 
     122                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra) 
    125123               ENDIF 
    126124            ENDDO 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90

    r11590 r11603  
    3030 
    3131#if ! defined PPR_LIB 
    32    SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
     32   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, kjpk_in, kjpk_out, kn_var)       
    3333      !!---------------------------------------------------------------------- 
    3434      !!                    ***  ROUTINE  reconstructandremap *** 
    … …  
    4747      !!                Give references if exist otherwise suppress these lines 
    4848      !!----------------------------------------------------------------------- 
    49       INTEGER N, Nout 
    50       REAL(wp) tabin(N), tabout(Nout) 
    51       REAL(wp) hin(N), hout(Nout) 
    52       REAL(wp) coeffremap(N,3),zwork(N,3) 
    53       REAL(wp) zwork2(N+1,3) 
    54       INTEGER jk 
    55       DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8   
    56       REAL(wp) q,q01,q02,q001,q002,q0 
    57       REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 
    58       REAL(wp),PARAMETER :: dpthin = 1.D-3 
    59       INTEGER :: k1, kbox, ktop, ka, kbot 
    60       REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
     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) 
    6163      !!----------------------------------------------------------------------- 
    6264 
    63       z_win(1)=0.; z_wout(1)= 0. 
    64       DO jk=1,N 
    65          z_win(jk+1)=z_win(jk)+hin(jk) 
    66       ENDDO  
     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  
    6769       
    68       DO jk=1,Nout 
    69          z_wout(jk+1)=z_wout(jk)+hout(jk)        
    70       ENDDO        
    71  
    72       DO jk=2,N 
    73          zwork(jk,1)=1./(hin(jk-1)+hin(jk)) 
    74       ENDDO 
    75          
    76       DO jk=2,N-1 
    77          q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) 
    78          zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) 
    79          zwork(jk,3)=q0 
    80       ENDDO        
     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./(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        
    8183      
    82       DO jk= 2,N 
    83          zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) 
    84       ENDDO 
    85          
    86       coeffremap(:,1) = tabin(:) 
     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) 
    8791  
    88       DO jk=2,N-1 
    89          q001 = hin(jk)*zwork2(jk+1,1) 
    90          q002 = hin(jk)*zwork2(jk,1)         
    91          IF (q001*q002 < 0) then 
    92             q001 = 0. 
    93             q002 = 0. 
    94          ENDIF 
    95          q=zwork(jk,2) 
    96          q01=q*zwork2(jk+1,1) 
    97          q02=q*zwork2(jk,1) 
    98          IF (abs(q001) > abs(q02)) q001 = q02 
    99          IF (abs(q002) > abs(q01)) q002 = q01 
    100          
    101          q=(q001-q002)*zwork(jk,3) 
    102          q001=q001-q*hin(jk+1) 
    103          q002=q002+q*hin(jk-1) 
    104          
    105          coeffremap(jk,3)=coeffremap(jk,1)+q001 
    106          coeffremap(jk,2)=coeffremap(jk,1)-q002 
    107          
    108          zwork2(jk,1)=(2.*q001-q002)**2 
    109          zwork2(jk,2)=(2.*q002-q001)**2 
    110       ENDDO 
    111          
    112       DO jk=1,N 
    113          IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 
    114             coeffremap(jk,3) = coeffremap(jk,1) 
    115             coeffremap(jk,2) = coeffremap(jk,1) 
    116             zwork2(jk,1) = 0. 
    117             zwork2(jk,2) = 0. 
    118          ENDIF 
    119       ENDDO 
    120          
    121       DO jk=2,N 
    122          q002=max(zwork2(jk-1,2),dsmll) 
    123          q001=max(zwork2(jk,1),dsmll) 
    124          zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) 
    125       ENDDO 
    126          
    127       zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
    128       zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
     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) 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*coeffremap(1,1)-zwork2(2,3) 
     132         zwork2(kjpk_in+1,3)=2*coeffremap(kjpk_in,1)-zwork2(kjpk_in,3) 
    129133  
    130       DO jk=1,N 
    131          q01=zwork2(jk+1,3)-coeffremap(jk,1) 
    132          q02=coeffremap(jk,1)-zwork2(jk,3) 
    133          q001=2.*q01 
    134          q002=2.*q02 
    135          IF (q01*q02<0) then 
    136             q01=0. 
    137             q02=0. 
    138          ELSEIF (abs(q01)>abs(q002)) then 
    139             q01=q002 
    140          ELSEIF (abs(q02)>abs(q001)) then 
    141             q02=q001 
    142          ENDIF 
    143          coeffremap(jk,2)=coeffremap(jk,1)-q02 
    144          coeffremap(jk,3)=coeffremap(jk,1)+q01 
    145       ENDDO 
    146  
    147       zbot=0.0 
    148       kbot=1 
    149       DO jk=1,Nout 
    150          ztop=zbot  !top is bottom of previous layer 
    151          ktop=kbot 
    152          IF     (ztop.GE.z_win(ktop+1)) then 
    153             ktop=ktop+1 
    154          ENDIF 
    155          
    156          zbot=z_wout(jk+1) 
    157          zthk=zbot-ztop 
    158  
    159          IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 
    160  
    161             kbot=ktop 
    162             DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
    163                kbot=kbot+1 
    164             ENDDO 
    165             zbox=zbot 
    166             DO k1= jk+1,Nout 
    167                IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN 
    168                   exit !thick layer 
     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.*q01 
     138            q002=2.*q02 
     139            IF (q01*q02<0) 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.-2.*(q001+q002) 
     193                        q02=q01-1.+(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  
    169204               ELSE 
    170                   zbox=z_wout(k1+1)  !include thin adjacent layers 
    171                   IF(zbox.EQ.z_wout(Nout+1)) THEN 
    172                      exit !at bottom 
     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 
    173227                  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) 
    174256               ENDIF 
    175             ENDDO 
    176             zthk=zbox-ztop 
    177  
    178             kbox=ktop 
    179             DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
    180                kbox=kbox+1 
    181             ENDDO 
    182            
    183             IF(ktop.EQ.kbox) THEN 
    184                IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 
    185                   IF(hin(kbox).GT.dpthin) THEN 
    186                      q001 = (zbox-z_win(kbox))/hin(kbox) 
    187                      q002 = (ztop-z_win(kbox))/hin(kbox) 
    188                      q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
    189                      q02=q01-1.+(q001+q002) 
    190                      q0=1.-q01-q02 
    191                   ELSE 
    192                      q0 = 1.0 
    193                      q01 = 0. 
    194                      q02 = 0. 
    195                   ENDIF 
    196                   tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
    197                ELSE 
    198                   tabout(jk) = tabin(kbox) 
    199                ENDIF  
    200             ELSE 
    201                IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 
    202                   ka = jk 
    203                ELSEIF (kbox-ktop.GE.3) THEN 
    204                   ka = (kbox+ktop)/2 
    205                ELSEIF (hin(ktop).GE.hin(kbox)) THEN 
    206                   ka = ktop 
    207                ELSE 
    208                   ka = kbox 
    209                ENDIF !choose ka 
    210      
    211                offset=coeffremap(ka,1) 
    212      
    213                qtop = z_win(ktop+1)-ztop !partial layer thickness 
    214                IF(hin(ktop).GT.dpthin) THEN 
    215                   q=(ztop-z_win(ktop))/hin(ktop) 
    216                   q01=q*(q-1.) 
    217                   q02=q01+q 
    218                   q0=1-q01-q02             
    219                ELSE 
    220                   q0 = 1. 
    221                   q01 = 0. 
    222                   q02 = 0. 
    223                ENDIF 
    224                 
    225                tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
    226      
    227                DO k1= ktop+1,kbox-1 
    228                   tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
    229                ENDDO !k1 
    230                 
    231                qbot = zbox-z_win(kbox) !partial layer thickness 
    232                IF(hin(kbox).GT.dpthin) THEN 
    233                   q=qbot/hin(kbox) 
    234                   q01=(q-1.)**2 
    235                   q02=q01-1.+q 
    236                   q0=1-q01-q02                             
    237                ELSE 
    238                   q0 = 1.0 
    239                   q01 = 0. 
    240                   q02 = 0. 
    241                ENDIF 
    242                
    243                tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
    244                 
    245                rpsum=1.0d0/zthk 
    246                tabout(jk)=offset+tsum*rpsum 
    247                   
    248             ENDIF !single or multiple layers 
    249          ELSE 
    250             IF (jk==1) THEN 
    251                write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 
    252             ENDIF 
    253             tabout(jk) = tabout(jk-1) 
     257               ptout(jk,jn) = ptout(jk-1,jn) 
    254258              
    255          ENDIF !normal:thin layer 
    256       ENDDO !jk 
     259            ENDIF !normal:thin layer 
     260         ENDDO !jk 
     261 
     262      END DO ! loop over variables 
    257263             
    258264   END SUBROUTINE reconstructandremap 
    … …  
    260266#else 
    261267 
    262    SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, jkin, jkout) 
     268   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, kjpk_in, kjpk_out, kn_var) 
    263269      !!---------------------------------------------------------------------- 
    264       !!                    ***  ROUTINE  reconstructandremap *** 
    265       !! 
    266       !! ** Purpose :   Brief description of the routine 
    267       !! 
    268       !! ** Method  :   description of the methodoloy used to achieve the 
    269       !!                objectives of the routine. Be as clear as possible! 
    270       !! 
    271       !! ** Action  : - first action (share memory array/varible modified 
    272       !!                in this routine 
    273       !!              - second action ..... 
    274       !!              - ..... 
    275       !! 
    276       !! References :   Author et al., Short_name_review, Year 
    277       !!                Give references if exist otherwise suppress these lines 
     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 
    278282      !!----------------------------------------------------------------------- 
    279       INTEGER , INTENT(in   )                     ::   jkin, jkout       !  
    280       REAL(wp), INTENT(in   ), DIMENSION(jkin)    ::   phin, ptin        ! 
    281       REAL(wp), INTENT(in   ), DIMENSION(jkout)   ::   phout             ! 
    282       REAL(wp), INTENT(inout), DIMENSION(jkout)   ::   ptout             ! 
     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 
    283290      ! 
    284       INTEGER, PARAMETER :: ndof = 1, nvar = 1 
    285       INTEGER  :: jk 
    286       REAL(wp) ::  zwin(jkin+1),   ztin(ndof, nvar, jkin) 
    287       REAL(wp) :: zwout(jkout+1), ztout(ndof, nvar, jkout) 
     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) 
    288295      TYPE(rmap_work) :: work 
    289296      TYPE(rmap_opts) :: opts 
    290       TYPE(rcon_ends) :: bc_l(nvar) 
    291       TYPE(rcon_ends) :: bc_r(nvar) 
     297      TYPE(rcon_ends) :: bc_l(kn_var) 
     298      TYPE(rcon_ends) :: bc_r(kn_var) 
    292299      !!-------------------------------------------------------------------- 
    293300      
    294301      ! Set interfaces and input data: 
    295302      zwin(1) = 0._wp 
    296       DO jk = 2, jkin + 1 
     303      DO jk = 2, kjpk_in + 1 
    297304         zwin(jk) = zwin(jk-1) + phin(jk-1)  
    298          ztin(ndof, nvar,jk-1) =  ptin(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 
    299311      END DO 
    300312 
    301313      zwout(1) = 0._wp 
    302       DO jk = 2, jkout + 1 
     314      DO jk = 2, kjpk_out + 1 
    303315         zwout(jk) = zwout(jk-1) + phout(jk-1)  
    304316      END DO 
    305317 
    306318      ! specify methods 
    307 !      opts%edge_meth = p3e_method     ! 3rd-order edge interp. 
    308 !      opts%cell_meth = ppm_method     ! PPM method in cells 
    309 !      opts%cell_lims = null_limit     ! no lim.      
    310       opts%edge_meth = p5e_method     ! 5th-order edge interp. 
    311       opts%cell_meth = pqm_method     ! PPM method in cells 
    312       opts%cell_lims = mono_limit     ! monotone limiter    
     319      opts%edge_meth = p3e_method     ! 3rd-order edge interp. 
     320      opts%cell_meth = ppm_method     ! PPM method in cells 
     321      opts%cell_lims = null_limit     ! no lim.      
     322!      opts%edge_meth = p5e_method     ! 5th-order edge interp. 
     323!      opts%cell_meth = pqm_method     ! PQM method in cells 
     324!      opts%cell_lims = mono_limit     ! monotone limiter    
    313325  
    314326      ! set boundary conditions 
    … …  
    319331 
    320332      ! init. method workspace 
    321       CALL work%init(jkin+1, nvar, opts) 
     333      CALL work%init(kjpk_in+1, kn_var, opts) 
    322334 
    323335      ! remap 
    324       CALL rmap1d(jkin+1, jkout+1, nvar, ndof, & 
    325       &           zwin, zwout, ztin, ztout,    & 
     336      CALL rmap1d(kjpk_in+1, kjpk_out+1, kn_var, ndof, & 
     337      &           zwin, zwout, ztin, ztout,            & 
    326338      &           bc_l, bc_r, work, opts) 
    327339 
    … …  
    329341      CALL work%free() 
    330342 
    331       DO jk =1, jkout 
    332          ptout(jk) = ztout(1,1,jk) 
     343      DO jn = 1, kn_var  
     344         DO jk = 1, kjpk_out 
     345            ptout(jk, jn) = ztout(1, jn, jk) 
     346         END DO 
    333347      END DO 
    334348             
Note: See TracChangeset for help on using the changeset viewer.