Ignore:
Timestamp:
2019-09-23T18:25:29+02:00 (14 months ago)
Author:
jchanut
Message:

#2222: 1) create remapping module (vremap) and integration of D. Engwirda piecewise polynomial recontruction package (PPR_LIB cpp key). 2) Various bug corrections with key_vertical activated.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90

    r11574 r11590  
    1717 
    1818   PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 
    19 #if defined key_vertical 
    20    PUBLIC reconstructandremap ! remapping routine 
    21 #endif    
     19   
    2220   !                                              !!* Namelist namagrif: AGRIF parameters 
    2321   LOGICAL , PUBLIC ::   ln_agrif_2way = .TRUE.    !: activate two way nesting  
     
    9593   END FUNCTION agrif_oce_alloc 
    9694 
    97 #if defined key_vertical 
    98    SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
    99       !!---------------------------------------------------------------------- 
    100       !!                ***  FUNCTION reconstructandremap  *** 
    101       !!---------------------------------------------------------------------- 
    102       IMPLICIT NONE 
    103       INTEGER N, Nout 
    104       REAL(wp) tabin(N), tabout(Nout) 
    105       REAL(wp) hin(N), hout(Nout) 
    106       REAL(wp) coeffremap(N,3),zwork(N,3) 
    107       REAL(wp) zwork2(N+1,3) 
    108       INTEGER jk 
    109       DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8   
    110       REAL(wp) q,q01,q02,q001,q002,q0 
    111       REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 
    112       REAL(wp),PARAMETER :: dpthin = 1.D-3 
    113       INTEGER :: k1, kbox, ktop, ka, kbot 
    114       REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
    115  
    116       z_win(1)=0.; z_wout(1)= 0. 
    117       DO jk=1,N 
    118          z_win(jk+1)=z_win(jk)+hin(jk) 
    119       ENDDO  
    120        
    121       DO jk=1,Nout 
    122          z_wout(jk+1)=z_wout(jk)+hout(jk)        
    123       ENDDO        
    124  
    125       DO jk=2,N 
    126          zwork(jk,1)=1./(hin(jk-1)+hin(jk)) 
    127       ENDDO 
    128          
    129       DO jk=2,N-1 
    130          q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) 
    131          zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) 
    132          zwork(jk,3)=q0 
    133       ENDDO        
    134       
    135       DO jk= 2,N 
    136          zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) 
    137       ENDDO 
    138          
    139       coeffremap(:,1) = tabin(:) 
    140   
    141       DO jk=2,N-1 
    142          q001 = hin(jk)*zwork2(jk+1,1) 
    143          q002 = hin(jk)*zwork2(jk,1)         
    144          IF (q001*q002 < 0) then 
    145             q001 = 0. 
    146             q002 = 0. 
    147          ENDIF 
    148          q=zwork(jk,2) 
    149          q01=q*zwork2(jk+1,1) 
    150          q02=q*zwork2(jk,1) 
    151          IF (abs(q001) > abs(q02)) q001 = q02 
    152          IF (abs(q002) > abs(q01)) q002 = q01 
    153          
    154          q=(q001-q002)*zwork(jk,3) 
    155          q001=q001-q*hin(jk+1) 
    156          q002=q002+q*hin(jk-1) 
    157          
    158          coeffremap(jk,3)=coeffremap(jk,1)+q001 
    159          coeffremap(jk,2)=coeffremap(jk,1)-q002 
    160          
    161          zwork2(jk,1)=(2.*q001-q002)**2 
    162          zwork2(jk,2)=(2.*q002-q001)**2 
    163       ENDDO 
    164          
    165       DO jk=1,N 
    166          IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 
    167             coeffremap(jk,3) = coeffremap(jk,1) 
    168             coeffremap(jk,2) = coeffremap(jk,1) 
    169             zwork2(jk,1) = 0. 
    170             zwork2(jk,2) = 0. 
    171          ENDIF 
    172       ENDDO 
    173          
    174       DO jk=2,N 
    175          q002=max(zwork2(jk-1,2),dsmll) 
    176          q001=max(zwork2(jk,1),dsmll) 
    177          zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) 
    178       ENDDO 
    179          
    180       zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
    181       zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
    182   
    183       DO jk=1,N 
    184          q01=zwork2(jk+1,3)-coeffremap(jk,1) 
    185          q02=coeffremap(jk,1)-zwork2(jk,3) 
    186          q001=2.*q01 
    187          q002=2.*q02 
    188          IF (q01*q02<0) then 
    189             q01=0. 
    190             q02=0. 
    191          ELSEIF (abs(q01)>abs(q002)) then 
    192             q01=q002 
    193          ELSEIF (abs(q02)>abs(q001)) then 
    194             q02=q001 
    195          ENDIF 
    196          coeffremap(jk,2)=coeffremap(jk,1)-q02 
    197          coeffremap(jk,3)=coeffremap(jk,1)+q01 
    198       ENDDO 
    199  
    200       zbot=0.0 
    201       kbot=1 
    202       DO jk=1,Nout 
    203          ztop=zbot  !top is bottom of previous layer 
    204          ktop=kbot 
    205          IF     (ztop.GE.z_win(ktop+1)) then 
    206             ktop=ktop+1 
    207          ENDIF 
    208          
    209          zbot=z_wout(jk+1) 
    210          zthk=zbot-ztop 
    211  
    212          IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 
    213  
    214             kbot=ktop 
    215             DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
    216                kbot=kbot+1 
    217             ENDDO 
    218             zbox=zbot 
    219             DO k1= jk+1,Nout 
    220                IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN 
    221                   exit !thick layer 
    222                ELSE 
    223                   zbox=z_wout(k1+1)  !include thin adjacent layers 
    224                   IF(zbox.EQ.z_wout(Nout+1)) THEN 
    225                      exit !at bottom 
    226                   ENDIF 
    227                ENDIF 
    228             ENDDO 
    229             zthk=zbox-ztop 
    230  
    231             kbox=ktop 
    232             DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
    233                kbox=kbox+1 
    234             ENDDO 
    235            
    236             IF(ktop.EQ.kbox) THEN 
    237                IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 
    238                   IF(hin(kbox).GT.dpthin) THEN 
    239                      q001 = (zbox-z_win(kbox))/hin(kbox) 
    240                      q002 = (ztop-z_win(kbox))/hin(kbox) 
    241                      q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
    242                      q02=q01-1.+(q001+q002) 
    243                      q0=1.-q01-q02 
    244                   ELSE 
    245                      q0 = 1.0 
    246                      q01 = 0. 
    247                      q02 = 0. 
    248                   ENDIF 
    249                   tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
    250                ELSE 
    251                   tabout(jk) = tabin(kbox) 
    252                ENDIF  
    253             ELSE 
    254                IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 
    255                   ka = jk 
    256                ELSEIF (kbox-ktop.GE.3) THEN 
    257                   ka = (kbox+ktop)/2 
    258                ELSEIF (hin(ktop).GE.hin(kbox)) THEN 
    259                   ka = ktop 
    260                ELSE 
    261                   ka = kbox 
    262                ENDIF !choose ka 
    263      
    264                offset=coeffremap(ka,1) 
    265      
    266                qtop = z_win(ktop+1)-ztop !partial layer thickness 
    267                IF(hin(ktop).GT.dpthin) THEN 
    268                   q=(ztop-z_win(ktop))/hin(ktop) 
    269                   q01=q*(q-1.) 
    270                   q02=q01+q 
    271                   q0=1-q01-q02             
    272                ELSE 
    273                   q0 = 1. 
    274                   q01 = 0. 
    275                   q02 = 0. 
    276                ENDIF 
    277                 
    278                tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
    279      
    280                DO k1= ktop+1,kbox-1 
    281                   tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
    282                ENDDO !k1 
    283                 
    284                qbot = zbox-z_win(kbox) !partial layer thickness 
    285                IF(hin(kbox).GT.dpthin) THEN 
    286                   q=qbot/hin(kbox) 
    287                   q01=(q-1.)**2 
    288                   q02=q01-1.+q 
    289                   q0=1-q01-q02                             
    290                ELSE 
    291                   q0 = 1.0 
    292                   q01 = 0. 
    293                   q02 = 0. 
    294                ENDIF 
    295                
    296                tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
    297                 
    298                rpsum=1.0d0/zthk 
    299                tabout(jk)=offset+tsum*rpsum 
    300                   
    301             ENDIF !single or multiple layers 
    302          ELSE 
    303             IF (jk==1) THEN 
    304                write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 
    305             ENDIF 
    306             tabout(jk) = tabout(jk-1) 
    307               
    308          ENDIF !normal:thin layer 
    309       ENDDO !jk 
    310              
    311       return 
    312       end subroutine reconstructandremap 
    313 #endif 
    314  
    31595#endif 
    31696   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.