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 8855 – NEMO

Changeset 8855


Ignore:
Timestamp:
2017-11-30T12:53:57+01:00 (7 years ago)
Author:
timgraham
Message:

Tidying of reconstructandremap subroutine

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90

    r8135 r8855  
    105105   END FUNCTION agrif_oce_alloc 
    106106 
    107       subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
    108       implicit none 
    109       integer N, Nout 
    110       real tabin(N), tabout(Nout) 
    111       real hin(N), hout(Nout) 
    112       real coeffremap(N,3),zwork(N,3) 
    113       real zwork2(N+1,3) 
    114       integer k 
    115       double precision, parameter :: dsmll=1.0d-8   
    116       real q,q01,q02,q001,q002,q0 
    117       real z_win(1:N+1), z_wout(1:Nout+1) 
    118       real,parameter :: dpthin = 1.D-3 
    119       integer :: k1, kbox, ktop, ka, kbot 
    120       real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
     107   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
     108      !!---------------------------------------------------------------------- 
     109      !!                ***  FUNCTION reconstructandremap  *** 
     110      !!---------------------------------------------------------------------- 
     111      IMPLICIT NONE 
     112      INTEGER N, Nout 
     113      REAL(wp) tabin(N), tabout(Nout) 
     114      REAL(wp) hin(N), hout(Nout) 
     115      REAL(wp) coeffremap(N,3),zwork(N,3) 
     116      REAL(wp) zwork2(N+1,3) 
     117      INTEGER jk 
     118      DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8   
     119      REAL(wp) q,q01,q02,q001,q002,q0 
     120      REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 
     121      REAL(wp),PARAMETER :: dpthin = 1.D-3 
     122      INTEGER :: k1, kbox, ktop, ka, kbot 
     123      REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
    121124 
    122125      z_win(1)=0.; z_wout(1)= 0. 
    123       do k=1,N 
    124        z_win(k+1)=z_win(k)+hin(k) 
    125       enddo  
     126      DO jk=1,N 
     127         z_win(jk+1)=z_win(jk)+hin(jk) 
     128      ENDDO  
    126129       
    127       do k=1,Nout 
    128        z_wout(k+1)=z_wout(k)+hout(k)        
    129       enddo        
    130  
    131         do k=2,N 
    132           zwork(k,1)=1./(hin(k-1)+hin(k)) 
    133         enddo 
    134          
    135         do k=2,N-1 
    136           q0 = 1./(hin(k-1)+hin(k)+hin(k+1)) 
    137           zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1) 
    138           zwork(k,3)=q0 
    139         enddo        
    140        
    141         do k= 2,N 
    142         zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1)) 
    143         enddo 
    144          
    145         coeffremap(:,1) = tabin(:) 
     130      DO jk=1,Nout 
     131         z_wout(jk+1)=z_wout(jk)+hout(jk)        
     132      ENDDO        
     133 
     134      DO jk=2,N 
     135         zwork(jk,1)=1./(hin(jk-1)+hin(jk)) 
     136      ENDDO 
     137         
     138      DO jk=2,N-1 
     139         q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) 
     140         zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) 
     141         zwork(jk,3)=q0 
     142      ENDDO        
     143      
     144      DO jk= 2,N 
     145         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) 
     146      ENDDO 
     147         
     148      coeffremap(:,1) = tabin(:) 
    146149  
    147          do k=2,N-1 
    148         q001 = hin(k)*zwork2(k+1,1) 
    149         q002 = hin(k)*zwork2(k,1)         
    150         if (q001*q002 < 0) then 
    151           q001 = 0. 
    152           q002 = 0. 
    153         endif 
    154         q=zwork(k,2) 
    155         q01=q*zwork2(k+1,1) 
    156         q02=q*zwork2(k,1) 
    157         if (abs(q001) > abs(q02)) q001 = q02 
    158         if (abs(q002) > abs(q01)) q002 = q01 
    159          
    160         q=(q001-q002)*zwork(k,3) 
    161         q001=q001-q*hin(k+1) 
    162         q002=q002+q*hin(k-1) 
    163          
    164         coeffremap(k,3)=coeffremap(k,1)+q001 
    165         coeffremap(k,2)=coeffremap(k,1)-q002 
    166          
    167         zwork2(k,1)=(2.*q001-q002)**2 
    168         zwork2(k,2)=(2.*q002-q001)**2 
    169         enddo 
    170          
    171         do k=1,N 
    172         if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then 
    173         coeffremap(k,3) = coeffremap(k,1) 
    174         coeffremap(k,2) = coeffremap(k,1) 
    175         zwork2(k,1) = 0. 
    176         zwork2(k,2) = 0. 
    177         endif 
    178         enddo 
    179          
    180         do k=2,N 
    181         q002=max(zwork2(k-1,2),dsmll) 
    182         q001=max(zwork2(k,1),dsmll) 
    183         zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002) 
    184         enddo 
    185          
    186         zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
    187         zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
     150      DO jk=2,N-1 
     151         q001 = hin(jk)*zwork2(jk+1,1) 
     152         q002 = hin(jk)*zwork2(jk,1)         
     153         IF (q001*q002 < 0) then 
     154            q001 = 0. 
     155            q002 = 0. 
     156         ENDIF 
     157         q=zwork(jk,2) 
     158         q01=q*zwork2(jk+1,1) 
     159         q02=q*zwork2(jk,1) 
     160         IF (abs(q001) > abs(q02)) q001 = q02 
     161         IF (abs(q002) > abs(q01)) q002 = q01 
     162         
     163         q=(q001-q002)*zwork(jk,3) 
     164         q001=q001-q*hin(jk+1) 
     165         q002=q002+q*hin(jk-1) 
     166         
     167         coeffremap(jk,3)=coeffremap(jk,1)+q001 
     168         coeffremap(jk,2)=coeffremap(jk,1)-q002 
     169         
     170         zwork2(jk,1)=(2.*q001-q002)**2 
     171         zwork2(jk,2)=(2.*q002-q001)**2 
     172      ENDDO 
     173         
     174      DO jk=1,N 
     175         IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 
     176            coeffremap(jk,3) = coeffremap(jk,1) 
     177            coeffremap(jk,2) = coeffremap(jk,1) 
     178            zwork2(jk,1) = 0. 
     179            zwork2(jk,2) = 0. 
     180         ENDIF 
     181      ENDDO 
     182         
     183      DO jk=2,N 
     184         q002=max(zwork2(jk-1,2),dsmll) 
     185         q001=max(zwork2(jk,1),dsmll) 
     186         zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) 
     187      ENDDO 
     188         
     189      zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
     190      zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
    188191  
    189         do k=1,N 
    190         q01=zwork2(k+1,3)-coeffremap(k,1) 
    191         q02=coeffremap(k,1)-zwork2(k,3) 
    192         q001=2.*q01 
    193         q002=2.*q02 
    194         if (q01*q02<0) then 
    195           q01=0. 
    196           q02=0. 
    197         elseif (abs(q01)>abs(q002)) then 
    198           q01=q002 
    199         elseif (abs(q02)>abs(q001)) then 
    200           q02=q001 
    201         endif 
    202         coeffremap(k,2)=coeffremap(k,1)-q02 
    203         coeffremap(k,3)=coeffremap(k,1)+q01 
    204         enddo 
     192      DO jk=1,N 
     193         q01=zwork2(jk+1,3)-coeffremap(jk,1) 
     194         q02=coeffremap(jk,1)-zwork2(jk,3) 
     195         q001=2.*q01 
     196         q002=2.*q02 
     197         IF (q01*q02<0) then 
     198            q01=0. 
     199            q02=0. 
     200         ELSEIF (abs(q01)>abs(q002)) then 
     201            q01=q002 
     202         ELSEIF (abs(q02)>abs(q001)) then 
     203            q02=q001 
     204         ENDIF 
     205         coeffremap(jk,2)=coeffremap(jk,1)-q02 
     206         coeffremap(jk,3)=coeffremap(jk,1)+q01 
     207      ENDDO 
    205208 
    206209      zbot=0.0 
    207210      kbot=1 
    208       do k=1,Nout 
    209         ztop=zbot  !top is bottom of previous layer 
    210         ktop=kbot 
    211         if     (ztop.ge.z_win(ktop+1)) then 
    212           ktop=ktop+1 
    213         endif 
    214          
    215         zbot=z_wout(k+1) 
    216         zthk=zbot-ztop 
    217  
    218         if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then 
    219  
    220           kbot=ktop 
    221           do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
    222             kbot=kbot+1 
    223           enddo 
    224           zbox=zbot 
    225           do k1= k+1,Nout 
    226             if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then 
    227               exit !thick layer 
    228             else 
    229               zbox=z_wout(k1+1)  !include thin adjacent layers 
    230               if     (zbox.eq.z_wout(Nout+1)) then 
    231                 exit !at bottom 
    232               endif 
    233             endif 
    234           enddo 
    235           zthk=zbox-ztop 
    236  
    237           kbox=ktop 
    238           do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
    239             kbox=kbox+1 
    240           enddo 
     211      DO jk=1,Nout 
     212         ztop=zbot  !top is bottom of previous layer 
     213         ktop=kbot 
     214         IF     (ztop.GE.z_win(ktop+1)) then 
     215            ktop=ktop+1 
     216         ENDIF 
     217         
     218         zbot=z_wout(jk+1) 
     219         zthk=zbot-ztop 
     220 
     221         IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 
     222 
     223            kbot=ktop 
     224            DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
     225               kbot=kbot+1 
     226            ENDDO 
     227            zbox=zbot 
     228            DO k1= jk+1,Nout 
     229               IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN 
     230                  exit !thick layer 
     231               ELSE 
     232                  zbox=z_wout(k1+1)  !include thin adjacent layers 
     233                  IF(zbox.EQ.z_wout(Nout+1)) THEN 
     234                     exit !at bottom 
     235                  ENDIF 
     236               ENDIF 
     237            ENDDO 
     238            zthk=zbox-ztop 
     239 
     240            kbox=ktop 
     241            DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
     242               kbox=kbox+1 
     243            ENDDO 
    241244           
    242           if     (ktop.eq.kbox) then 
    243  
    244  
    245             if     (z_wout(k)  .ne.z_win(kbox)   .or.z_wout(k+1).ne.z_win(kbox+1)     ) then 
    246  
    247               if     (hin(kbox).gt.dpthin) then 
    248                 q001 = (zbox-z_win(kbox))/hin(kbox) 
    249                 q002 = (ztop-z_win(kbox))/hin(kbox) 
    250                 q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
    251                 q02=q01-1.+(q001+q002) 
    252                 q0=1.-q01-q02 
    253               else 
    254                 q0 = 1.0 
    255                 q01 = 0. 
    256                 q02 = 0. 
    257               endif 
    258           tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
     245            IF(ktop.EQ.kbox) THEN 
     246               IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 
     247                  IF(hin(kbox).GT.dpthin) THEN 
     248                     q001 = (zbox-z_win(kbox))/hin(kbox) 
     249                     q002 = (ztop-z_win(kbox))/hin(kbox) 
     250                     q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
     251                     q02=q01-1.+(q001+q002) 
     252                     q0=1.-q01-q02 
     253                  ELSE 
     254                     q0 = 1.0 
     255                     q01 = 0. 
     256                     q02 = 0. 
     257                  ENDIF 
     258                  tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
     259               ELSE 
     260                  tabout(jk) = tabin(kbox) 
     261               ENDIF  
     262            ELSE 
     263               IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 
     264                  ka = jk 
     265               ELSEIF (kbox-ktop.GE.3) THEN 
     266                  ka = (kbox+ktop)/2 
     267               ELSEIF (hin(ktop).GE.hin(kbox)) THEN 
     268                  ka = ktop 
     269               ELSE 
     270                  ka = kbox 
     271               ENDIF !choose ka 
     272     
     273               offset=coeffremap(ka,1) 
     274     
     275               qtop = z_win(ktop+1)-ztop !partial layer thickness 
     276               IF(hin(ktop).GT.dpthin) THEN 
     277                  q=(ztop-z_win(ktop))/hin(ktop) 
     278                  q01=q*(q-1.) 
     279                  q02=q01+q 
     280                  q0=1-q01-q02             
     281               ELSE 
     282                  q0 = 1. 
     283                  q01 = 0. 
     284                  q02 = 0. 
     285               ENDIF 
     286                
     287               tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
     288     
     289               DO k1= ktop+1,kbox-1 
     290                  tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
     291               ENDDO !k1 
     292                
     293               qbot = zbox-z_win(kbox) !partial layer thickness 
     294               IF(hin(kbox).GT.dpthin) THEN 
     295                  q=qbot/hin(kbox) 
     296                  q01=(q-1.)**2 
     297                  q02=q01-1.+q 
     298                  q0=1-q01-q02                             
     299               ELSE 
     300                  q0 = 1.0 
     301                  q01 = 0. 
     302                  q02 = 0. 
     303               ENDIF 
    259304               
    260             else 
    261             tabout(k) = tabin(kbox) 
    262                
    263             endif  
    264  
    265           else 
    266  
    267             if     (ktop.le.k .and. kbox.ge.k) then 
    268               ka = k 
    269             elseif (kbox-ktop.ge.3) then 
    270               ka = (kbox+ktop)/2 
    271             elseif (hin(ktop).ge.hin(kbox)) then 
    272               ka = ktop 
    273             else 
    274               ka = kbox 
    275             endif !choose ka 
    276  
    277             offset=coeffremap(ka,1) 
    278  
    279             qtop = z_win(ktop+1)-ztop !partial layer thickness 
    280             if     (hin(ktop).gt.dpthin) then 
    281               q=(ztop-z_win(ktop))/hin(ktop) 
    282               q01=q*(q-1.) 
    283               q02=q01+q 
    284               q0=1-q01-q02             
    285             else 
    286               q0 = 1. 
    287               q01 = 0. 
    288               q02 = 0. 
    289             endif 
    290              
    291             tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
    292  
    293             do k1= ktop+1,kbox-1 
    294               tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
    295             enddo !k1 
    296  
    297              
    298             qbot = zbox-z_win(kbox) !partial layer thickness 
    299             if     (hin(kbox).gt.dpthin) then 
    300               q=qbot/hin(kbox) 
    301               q01=(q-1.)**2 
    302               q02=q01-1.+q 
    303               q0=1-q01-q02                             
    304             else 
    305               q0 = 1.0 
    306               q01 = 0. 
    307               q02 = 0. 
    308             endif 
    309             
    310             tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
    311              
    312             rpsum=1.0d0/zthk 
    313               tabout(k)=offset+tsum*rpsum 
    314                
    315           endif !single or multiple layers 
    316         else 
    317         if (k==1) then 
    318         write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1) 
    319         endif 
    320          tabout(k) = tabout(k-1) 
    321            
    322         endif !normal:thin layer 
    323       enddo !k 
     305               tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
     306                
     307               rpsum=1.0d0/zthk 
     308               tabout(jk)=offset+tsum*rpsum 
     309                  
     310            ENDIF !single or multiple layers 
     311         ELSE 
     312            IF (jk==1) THEN 
     313               write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 
     314            ENDIF 
     315            tabout(jk) = tabout(jk-1) 
     316              
     317         ENDIF !normal:thin layer 
     318      ENDDO !jk 
    324319             
    325320      return 
Note: See TracChangeset for help on using the changeset viewer.