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 2908 for branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 – NEMO

Ignore:
Timestamp:
2011-10-13T07:52:23+02:00 (13 years ago)
Author:
cbricaud
Message:

add a s-coordinate case for interp subroutine

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    r2864 r2908  
    2121  !!---------------------------------------------------------------------- 
    2222  !!---------------------------------------------------------------------- 
    23   !!   dia_dct      :  compute the transport trough a sec. 
     23  !!   dia_dct      :  compute the transport through a sec. 
    2424  !!   dia_dct_init :  read namelist. 
    2525  !!   readsec      :  read sections description and pathway 
    … …  
    7878 
    7979  TYPE SECTION 
    80      CHARACTER(len=60)                            :: name           ! name of the sec 
    81      LOGICAL                                      :: llstrpond      ! true if you want the computation of salinity and 
    82                                                                     ! temperature balanced by the transport 
    83      LOGICAL                                      :: ll_ice_section ! icesurf and icevol computation 
    84      LOGICAL                                      :: ll_date_line   ! = T if the section crosses the date-line 
    85      TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec       ! longitude and latitude of the extremities of the sec 
    86      INTEGER                                      :: nb_class       ! number of boundaries for density classes 
    87      INTEGER, DIMENSION(nb_point_max)             :: direction      ! vector direction of the point in the section 
    88      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname      ! caracteristics of the class 
    89      REAL(wp), DIMENSION(nb_class_max)            :: zsigi        ,&! in-situ   density classes    (99 if you don't want) 
    90                                                      zsigp        ,&! potential density classes    (99 if you don't want) 
    91                                                      zsal         ,&! salinity classes   (99 if you don't want) 
    92                                                      ztem         ,&! temperature classes(99 if you don't want) 
    93                                                      zlay           ! level classes      (99 if you don't want) 
     80     CHARACTER(len=60)                            :: name              ! name of the sec 
     81     LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and 
     82                                                                       ! heat transports 
     83     LOGICAL                                      :: ll_ice_section    ! icesurf and icevol computation 
     84     LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line 
     85     TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec 
     86     INTEGER                                      :: nb_class          ! number of boundaries for density classes 
     87     INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section 
     88     CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! caracteristics of the class 
     89     REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want) 
     90                                                     zsigp           ,&! potential density classes    (99 if you don't want) 
     91                                                     zsal            ,&! salinity classes   (99 if you don't want) 
     92                                                     ztem            ,&! temperature classes(99 if you don't want) 
     93                                                     zlay              ! level classes      (99 if you don't want) 
    9494     REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
    95      REAL(wp)                                         :: slopeSection  ! coeff directeur de la section 
    96      INTEGER                                          :: nb_point              ! nb de points de la section 
    97      TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint             ! list of point in the section 
     95     REAL(wp)                                         :: slopeSection  ! section's slopesection  
     96     INTEGER                                          :: nb_point      ! section's number of points 
     97     TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! section's list point 
    9898  END TYPE SECTION 
    9999 
    … …  
    161161     LOGICAL             :: lldebug =.FALSE.  !debug a section   
    162162     CHARACTER(len=160)  :: clfileout         !fileout name 
     163 
     164      
     165     INTEGER , DIMENSION(1):: ish                                       ! tmp array for mpp_sum 
     166     INTEGER , DIMENSION(3):: ish2                                      !   " 
     167     REAL(wp), DIMENSION(nb_sec_max*nb_type_class*nb_class_max):: zwork !   "   
     168     REAL(wp), DIMENSION(nb_sec_max,nb_type_class,nb_class_max):: zsum  !   " 
     169 
    163170     !!---------------------------------------------------------------------     
    164171 
    … …  
    183190           !Compute transport through section   
    184191           CALL transport(secs(jsec),lldebug)  
     192 
     193        ENDDO 
    185194              
    186            IF( MOD(kt,nn_dctwri)==0 )THEN 
    187  
    188               IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: write at kt = ",kt          
     195        IF( MOD(kt,nn_dctwri)==0 )THEN 
     196 
     197           IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: write at kt = ",kt          
    189198   
    190               !Write the transport 
     199           !Sum on all procs  
     200           IF( lk_mpp )THEN 
     201              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max  
     202              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/) 
     203              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO 
     204              zwork(:)= RESHAPE(zsum(:,:,:), ish ) 
     205              CALL mpp_sum(zwork, ish(1)) 
     206              zsum(:,:,:)= RESHAPE(zwork,ish2) 
     207              DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO 
     208           ENDIF 
     209 
     210           !Write the transport 
     211           DO jsec=1,nb_sec 
     212 
    191213              IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
    192214             
    … …  
    194216              secs(jsec)%transport(:,:)=0.   
    195217 
    196            ENDIF  
    197         ENDDO 
     218           ENDDO 
     219 
     220        ENDIF  
     221 
    198222     ENDIF 
    199223 
    … …  
    310334              ENDDO                   
    311335           ENDIF 
    312            !?IF( jsec==nn_secdebug .OR. nn_secdebug==-1  )THEN 
    313            !  DO jpt=1,iptglo 
    314            !      WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt) 
    315            !   ENDDO                   
    316            !ENDIF 
    317336  
    318337           !Now each proc selects only points that are in its domain: 
    … …  
    335354                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates 
    336355                 secs(jsec)%direction(iptloc) = directemp(jpt)               ! store local direction 
    337                  !WRITE(narea+200,*)'         # I J : ',iiloc,ijloc 
    338356              ENDIF 
    339357 
    … …  
    356374 
    357375           !remove redundant points between processors 
    358            !look NEMO documentation for more explanations 
    359 !           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE. 
    360            lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1)  ) lldebug = .TRUE. 
     376           !------------------------------------------ 
     377           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE. 
    361378           IF( iptloc .NE. 0 )THEN 
    362379              CALL removepoints(secs(jsec),'I','top_list',lldebug) 
    … …  
    409426                                 !where points could be redundant 
    410427                isgn          ,& !way of course in listpoint 
     428                ipoint          ,& !way of course in listpoint 
    411429                istart,iend      !first and last points selected in listpoint 
    412430     INTEGER :: jpoint   =0      !loop on list points 
    … …  
    443461 
    444462     jpoint=iextr+isgn 
    445      DO WHILE ( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point .AND. & 
    446              icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest) 
    447              jpoint=jpoint+isgn 
     463     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point .AND. & 
     464        icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest ) 
     465        jpoint=jpoint+isgn 
    448466     ENDDO 
     467 
    449468     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point 
    450469     ELSE                        ; istart=1        ; iend=jpoint+1 
    … …  
    489508     !! * Local variables 
    490509     INTEGER             :: jk,jseg,jclass,   &!loop on level/segment/classes  
    491                             iptegrid,         &!=0 if section is constant along i/j  
    492510                            isgnu  , isgnv     ! 
    493511     INTEGER :: ii, ij ! local integer 
    … …  
    506524 
    507525     TYPE(POINT_SECTION) :: k 
    508      INTEGER,DIMENSION(1):: ish 
    509      INTEGER,DIMENSION(2):: ish2 
    510526     REAL(wp),DIMENSION(nb_type_class,nb_class_max)::zsum 
    511      REAL(wp),DIMENSION(nb_type_class*nb_class_max)::zwork  ! temporary vector for mpp_sum 
    512527     !!-------------------------------------------------------- 
    513528 
    … …  
    520535     zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp 
    521536     zice_vol_pos  = 0._wp ; zice_vol_neg  = 0._wp 
    522  
    523      !compute iptegrid: iptegrid=0 if section is constant along i or j 
    524      iptegrid = 0 
    525      jseg     = 1 
    526      DO WHILE( (jseg .LT. MAX(sec%nb_point-1,0)) .AND. iptegrid .NE. 1 ) 
    527         IF(sec%direction(jseg) .NE. sec%direction(1)) iptegrid = 1     
    528         jseg = jseg+1 
    529      ENDDO  
    530      IF(lk_mpp) CALL mpp_sum(iptegrid) 
    531537 
    532538     !---------------------------! 
    … …  
    623629                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
    624630                 !----------------------------------------------! 
    625            IF(ld_debug)write(narea+200,*)k ,jk 
    626631  
    627632                 IF ( (    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.    & 
    … …  
    663668                    ENDIF 
    664669#endif 
    665            !IF(ld_debug)write(narea+200,*)k ,jk,zTnorm              
    666670                    !COMPUTE TRANSPORT  
    667671                    !zTnorm=transport through one cell for one class 
    … …  
    771775     ENDIF   !end of nb_inmesh=0 case 
    772776 
    773      !------------------! 
    774      !SUM ON ALL PROCS  ! 
    775      !------------------! 
    776      IF( lk_mpp )THEN 
    777          ish(1)  =  nb_type_class*nb_class_max ;  ish2(1)=nb_type_class ; ish2(2)=nb_class_max 
    778          zwork(:)= RESHAPE(zsum(:,:), ish ) 
    779          CALL mpp_sum(zwork, ish(1)) 
    780          zsum(:,:)= RESHAPE(zwork,ish2) 
    781      ENDIF 
    782       
    783777     !-------------------------------| 
    784778     !FINISH COMPUTING TRANSPORTS    | 
    … …  
    838832     !! Method: 
    839833     !!        1. Write volume transports in "volume_transport" 
    840      !!           Unit: Sv : Surface * Velocity / 1.e6  
     834     !!           Unit: Sv : area * Velocity / 1.e6  
    841835     !!  
    842836     !!        2. Write heat transports in "heat_transport" 
    843      !!           Unit: Peta W : Surface * Velocity * T * rhau * Cp / 1.e15 
     837     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15 
    844838     !!  
    845839     !!        3. Write salt transports in "salt_transport" 
    846      !!           Unit: Mega m^3 / s : Surface * Velocity * S / 1.e6 
     840     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6 
    847841     !! 
    848842     !!-------------------------------------------------------------  
    849843     !!arguments 
    850      INTEGER, INTENT(IN)          :: kt 
    851      TYPE(SECTION), INTENT(INOUT) :: sec 
    852      INTEGER ,INTENT(IN)          :: jsec 
     844     INTEGER, INTENT(IN)          :: kt         ! time-step 
     845     TYPE(SECTION), INTENT(INOUT) :: sec        ! section to write    
     846     INTEGER ,INTENT(IN)          :: jsec       ! section number 
    853847 
    854848     !!local declarations 
    855849     REAL(wp) ,DIMENSION(nb_type_class):: zsumclass 
    856      INTEGER               :: jcl,ji        !Dummy loop 
    857      CHARACTER(len=2)      :: classe 
    858      REAL(wp)              :: zbnd1,zbnd2 
     850     INTEGER               :: jcl,ji            ! Dummy loop 
     851     CHARACTER(len=2)      :: classe            ! Classname  
     852     REAL(wp)              :: zbnd1,zbnd2       ! Class bounds 
     853     REAL(wp)              :: zslope            ! section's slope coeff 
    859854     !!-------------------------------------------------------------  
    860855       
    861856     zsumclass(:)=0._wp 
    862                 
     857     zslope = sec%slopeSection        
     858 
     859  
    863860     DO jcl=1,MAX(1,sec%nb_class-1) 
    864861 
    … …  
    900897        ENDIF 
    901898        !temperature classes transports 
    902         IF( ( sec%ztem(jcl+1) .NE. 99.     ) .AND. & 
    903             ( sec%ztem(jcl+1) .NE. 99.     )       ) THEN 
     899        IF( ( sec%ztem(jcl+1) .NE. 99._wp     ) .AND. & 
     900            ( sec%ztem(jcl+1) .NE. 99._wp     )       ) THEN 
    904901           classe = 'T       ' 
    905902           zbnd1 = sec%ztem(jcl) 
    … …  
    908905                   
    909906        !write volume transport per class 
    910         WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, & 
     907        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,zslope, & 
    911908                              jcl,classe,zbnd1,zbnd2,& 
    912909                              sec%transport(1,jcl),sec%transport(2,jcl), & 
    … …  
    916913 
    917914           !write heat transport per class: 
    918            WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection,  & 
     915           WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,zslope,  & 
    919916                              jcl,classe,zbnd1,zbnd2,& 
    920917                              sec%transport(7,jcl)*1000._wp*rcp/1.e15,sec%transport(8,jcl)*1000._wp*rcp/1.e15, & 
    921918                              ( sec%transport(7,jcl)+sec%transport(8,jcl) )*1000._wp*rcp/1.e15 
    922919           !write salt transport per class 
    923            WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection,  & 
     920           WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,zslope,  & 
    924921                              jcl,classe,zbnd1,zbnd2,& 
    925922                              sec%transport(9,jcl)*1000._wp/1.e9,sec%transport(10,jcl)*1000._wp/1.e9,& 
    … …  
    934931 
    935932     !write total volume transport 
    936      WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, & 
     933     WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,zslope, & 
    937934                           jcl,"total",zbnd1,zbnd2,& 
    938935                           zsumclass(1),zsumclass(2),zsumclass(1)+zsumclass(2) 
    … …  
    941938 
    942939        !write total heat transport 
    943         WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection, & 
     940        WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,zslope, & 
    944941                           jcl,"total",zbnd1,zbnd2,& 
    945                            zsumclass(7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,(zsumclass(7)+zsumclass(8) )* 1000._wp*rcp/1.e15 
     942                           zsumclass(7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,& 
     943                           (zsumclass(7)+zsumclass(8) )* 1000._wp*rcp/1.e15 
    946944        !write total salt transport 
    947         WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection, & 
     945        WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,zslope, & 
    948946                           jcl,"total",zbnd1,zbnd2,& 
    949                            zsumclass(9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,(zsumclass(9)+zsumclass(10))*1000._wp/1.e9 
     947                           zsumclass(9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,& 
     948                           (zsumclass(9)+zsumclass(10))*1000._wp/1.e9 
    950949     ENDIF 
    951950 
    … …  
    953952     IF ( sec%ll_ice_section) THEN 
    954953        !write total ice volume transport 
    955         WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,& 
     954        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,zslope,& 
    956955                              jcl,"ice_vol",zbnd1,zbnd2,& 
    957956                              sec%transport(9,1),sec%transport(10,1),& 
    958957                              sec%transport(9,1)+sec%transport(10,1) 
    959958        !write total ice surface transport 
    960         WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,& 
     959        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,zslope,& 
    961960                              jcl,"ice_surf",zbnd1,zbnd2,& 
    962961                              sec%transport(11,1),sec%transport(12,1), & 
    … …  
    972971  !!---------------------------------------------------------------------- 
    973972  !! 
    974   !!   Purpose: compute Temperature/Salinity/density on U-point or V-point 
     973  !!   Purpose: compute Temperature/Salinity/density at U-point or V-point 
     974  !!   -------- 
    975975  !! 
    976   !!    |    I       |    I+1     | 
    977   !!    |            |            | 
    978   !!  --------------------------------     1.Compute Tbis: 
    979   !!    |            |            |          interpolate T(I,J,K) et T(I,J,K+1) 
    980   !!    |            |            | 
    981   !!    |            |            |        2. 1.Compute Temp/Salt/Rho à U point:  
    982   !! K  |    T       |            |          interpolate Tbis et T(I+U,J,K+1)   
    983   !!    |            |            | 
    984   !!    |            |            | 
    985   !!    |            |            | 
    986   !!  -------------------------------- 
    987   !!    |            |            | 
    988   !!    |            |            | 
    989   !!    |            |            | 
    990   !!K+1 |    Tbis    U     T      | 
    991   !!    |            |            | 
    992   !!    |    T       |            | 
    993   !!    |            |------------| 
    994   !!    |            | partials   | 
    995   !!    |            |  steps     | 
    996   !!  -------------------------------- 
     976  !!   Method: 
     977  !!   ------ 
     978  !! 
     979  !!   ====> full step and partial step 
     980  !!  
     981  !! 
     982  !!    |    I          |    I+1           |    Z=Temperature/Salinity/density at U-poinT 
     983  !!    |               |                  | 
     984  !!  ----------------------------------------  1. Veritcale interpolation: compute zbis 
     985  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1) 
     986  !!    |               |                  |       zbis =  
     987  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ] 
     988  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]  
     989  !!    |               |                  |  
     990  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point 
     991  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K)   
     992  !!    |     .         |                  | 
     993  !!    |     .         |                  |       Z= ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)  
     994  !!    |     .         |                  | 
     995  !!  ------------------------------------------ 
     996  !!    |     .         |                  | 
     997  !!    |     .         |                  | 
     998  !!    |     .         |                  | 
     999  !!K   |    zbis...... U ..ptab(I+1,J,K)  | 
     1000  !!    |     .         |                  | 
     1001  !!    | ptab(I,J,K)   |                  | 
     1002  !!    |               |------------------| 
     1003  !!    |               | partials         | 
     1004  !!    |               |  steps           | 
     1005  !!  ------------------------------------------- 
     1006  !!    <----zet1------><----zet2---------> 
     1007  !! 
     1008  !! 
     1009  !!   ====>  s-coordinate 
     1010  !!      
     1011  !!    |                |                  |   1. Compute distance between T1 and U points 
     1012  !!    |                |                  |      Compute distance between T2 and U points 
     1013  !!    |                | ptab(I+1,J,K)    |  
     1014  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point  
     1015  !!    |                |      ^           |     
     1016  !!    |                |      | zdep2     |     
     1017  !!    |                |      |           |     
     1018  !!    |       ^        U      v           | 
     1019  !!    |       |        |                  | 
     1020  !!    |       | zdep1  |                  |     
     1021  !!    |       v        |                  | 
     1022  !!    |      T1        |                  | 
     1023  !!    | ptab(I,J,K)    |                  |  
     1024  !!    |                |                  |  
     1025  !!    |                |                  |  
     1026  !! 
     1027  !!    <----zet1--------><----zet2---------> 
    9971028  !! 
    9981029  !!---------------------------------------------------------------------- 
    9991030  !*arguments 
    1000   INTEGER, INTENT(IN)                          :: ki, kj, kk 
    1001   CHARACTER(len=1), INTENT(IN)                 :: cd_point 
    1002   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab 
    1003   REAL(wp)                                     :: interp 
     1031  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point 
     1032  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V) 
     1033  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk ) 
     1034  REAL(wp)                                     :: interp       ! interpolated variable  
    10041035 
    10051036  !*local declations 
    1006   INTEGER :: ii1, ij1, ii2, ij2 
    1007   REAL(wp):: ze3w, zfse3, zmax1, zmax2, zbis 
     1037  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer 
     1038  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real 
     1039  REAL(wp):: zet1, zet2                                        ! weight for interpolation  
     1040  REAL(wp):: zdep1,zdep2                                       ! differences of depth 
    10081041  !!---------------------------------------------------------------------- 
    10091042 
    … …  
    10111044     ii1 = ki    ; ij1 = kj  
    10121045     ii2 = ki+1  ; ij2 = kj  
     1046 
     1047     zet1=e1t(ii1,ij1) 
     1048     zet2=e1t(ii2,ij2) 
     1049 
    10131050  ELSE ! cd_point=='V'  
    10141051     ii1 = ki    ; ij1 = kj  
    10151052     ii2 = ki    ; ij2 = kj+1   
     1053 
     1054     zet1=e2t(ii1,ij1) 
     1055     zet2=e2t(ii2,ij2) 
     1056 
    10161057  ENDIF 
    10171058 
     1059  IF( ln_sco )THEN   ! s-coordinate case 
     1060 
     1061     zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2  
     1062     zdep1 = fsdept(ii1,ij1,kk) - zdepu 
     1063     zdep2 = fsdept(ii2,ij2,kk) - zdepu 
     1064 
     1065     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) ) 
     1066     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) ) 
     1067   
     1068     ! result 
     1069     interp = umask(ii1,ij1,kk) * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )    
     1070 
     1071 
     1072  ELSE       ! full step or partial step case  
     1073 
    10181074#if defined key_vvl 
    10191075 
    1020   ze3w  = fsve3w(ii2,ij2,kk) - fsve3w(ii1,ij1,kk) 
    1021  
    1022   zfse3 = fsve3w(ii1,ij1,kk) * ( 1 + sshn(ii2,ij2) * mut(ii2,ij2,kk) ) 
    1023   zmax1 = ( fse3w(ii2,ij2,kk) - zfse3 )/ fse3w(ii2,ij2,kk) 
     1076     ze3w  = fsve3w(ii2,ij2,kk) - fsve3w(ii1,ij1,kk) 
     1077 
     1078     zfse3 = fsve3w(ii1,ij1,kk) * ( 1 + sshn(ii2,ij2) * mut(ii2,ij2,kk) ) 
     1079     zwgt1 = ( fse3w(ii2,ij2,kk) - zfse3 )/ fse3w(ii2,ij2,kk) 
    10241080   
    1025   zfse3 = fsve3w(ii2,ij2,kk) * ( 1 + sshn(ii1,ij1) * mut(ii1,ij1,kk) ) 
    1026   zmax2 = ( fse3w(ii1,ij1,kk) - zfse3 )/ fse3w(ii1,ij1,kk) 
     1081     zfse3 = fsve3w(ii2,ij2,kk) * ( 1 + sshn(ii1,ij1) * mut(ii1,ij1,kk) ) 
     1082     zwgt2 = ( fse3w(ii1,ij1,kk) - zfse3 )/ fse3w(ii1,ij1,kk) 
    10271083 
    10281084#else 
    10291085 
    1030   ze3w  = fse3w(ii2,ij2,kk)-fse3w(ii1,ij1,kk)  
    1031   zmax1 =  ze3w / fse3w(ii2,ij2,kk) 
    1032   zmax2 = -ze3w / fse3w(ii1,ij1,kk) 
     1086     ze3t  = fse3t(ii2,ij2,kk)-fse3t(ii1,ij1,kk)  
     1087     zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk) 
     1088     zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk) 
    10331089 
    10341090#endif 
    10351091 
    1036   IF(kk .NE. 1)THEN 
    1037  
    1038      IF( ze3w >= 0. )THEN  
    1039         !zbis 
    1040         zbis = ptab(ii2,ij2,kk) + zmax1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )  
    1041         ! result 
    1042         interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + zbis ) 
     1092     IF(kk .NE. 1)THEN 
     1093 
     1094        IF( ze3t >= 0. )THEN  
     1095           !zbis 
     1096           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )  
     1097           ! result 
     1098            interp = umask(ii1,ij1,kk) * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 
     1099        ELSE 
     1100           !zbis 
     1101           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) ) 
     1102           ! result 
     1103           interp = umask(ii1,ij1,kk) * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     1104        ENDIF     
     1105 
    10431106     ELSE 
    1044         !zbis 
    1045         zbis = ptab(ii1,ij1,kk) + zmax2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) ) 
    1046         ! result 
    1047         interp = umask(ii1,ij1,kk) * 0.5*( zbis +  ptab(ii2,ij2,kk) ) 
    1048      ENDIF     
    1049  
    1050   ELSE 
    1051      interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + ptab(ii2,ij2,kk) ) 
     1107        interp = umask(ii1,ij1,kk) * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     1108     ENDIF 
     1109 
    10521110  ENDIF 
     1111 
    10531112 
    10541113  END FUNCTION interp 
Note: See TracChangeset for help on using the changeset viewer.