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 7351 for branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 – NEMO

Ignore:
Timestamp:
2016-11-28T17:04:10+01:00 (7 years ago)
Author:
emanuelaclementi
Message:

ticket #1805 step 3: /2016/dev_INGV_UKMO_2016 aligned to the trunk at revision 7161

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    • Property svn:executable deleted
    r5505 r7351  
    11MODULE diadct 
    2   !!===================================================================== 
    3   !!                       ***  MODULE  diadct  *** 
    4   !! Ocean diagnostics: Compute the transport trough a sec. 
    5   !!=============================================================== 
    6   !! History :  
    7   !! 
    8   !!         original  : 02/99 (Y Drillet) 
    9   !!         addition  : 10/01 (Y Drillet, R Bourdalle Badie) 
    10   !!                   : 10/05 (M Laborie) F90 
    11   !!         addition  : 04/07 (G Garric) Ice sections 
    12   !!         bugfix    : 04/07 (C Bricaud) test on sec%nb_point 
    13   !!                                      initialisation of ztransp1,ztransp2,... 
    14   !!         nemo_v_3_4: 09/2011 (C Bricaud) 
    15   !! 
    16   !! 
    17   !!---------------------------------------------------------------------- 
     2   !!====================================================================== 
     3   !!                       ***  MODULE  diadct  *** 
     4   !! Ocean diagnostics: Compute the transport trough a sec. 
     5   !!====================================================================== 
     6   !! History :  OPA  ! 02/1999 (Y Drillet)  original code 
     7   !!                 ! 10/2001 (Y Drillet, R Bourdalle Badie) 
     8   !!   NEMO     1.0  ! 10/2005 (M Laborie) F90 
     9   !!            3.0  ! 04/2007 (G Garric) Ice sections 
     10   !!             -   ! 04/2007 (C Bricaud) test on sec%nb_point, initialisation of ztransp1,ztransp2,... 
     11   !!            3.4  ! 09/2011 (C Bricaud) 
     12   !!---------------------------------------------------------------------- 
    1813#if defined key_diadct 
    19   !!---------------------------------------------------------------------- 
    20   !!   'key_diadct' : 
    21   !!---------------------------------------------------------------------- 
    22   !!---------------------------------------------------------------------- 
    23   !!   dia_dct      :  Compute the transport through a sec. 
    24   !!   dia_dct_init :  Read namelist. 
    25   !!   readsec      :  Read sections description and pathway 
    26   !!   removepoints :  Remove points which are common to 2 procs 
    27   !!   transport    :  Compute transport for each sections 
    28   !!   dia_dct_wri  :  Write tranports results in ascii files 
    29   !!   interp       :  Compute temperature/salinity/density at U-point or V-point 
    30   !!    
    31   !!---------------------------------------------------------------------- 
    32   !! * Modules used 
    33   USE oce             ! ocean dynamics and tracers 
    34   USE dom_oce         ! ocean space and time domain 
    35   USE phycst          ! physical constants 
    36   USE in_out_manager  ! I/O manager 
    37   USE daymod          ! calendar 
    38   USE dianam          ! build name of file 
    39   USE lib_mpp         ! distributed memory computing library 
     14   !!---------------------------------------------------------------------- 
     15   !!   'key_diadct' : 
     16   !!---------------------------------------------------------------------- 
     17   !!---------------------------------------------------------------------- 
     18   !!   dia_dct      :  Compute the transport through a sec. 
     19   !!   dia_dct_init :  Read namelist. 
     20   !!   readsec      :  Read sections description and pathway 
     21   !!   removepoints :  Remove points which are common to 2 procs 
     22   !!   transport    :  Compute transport for each sections 
     23   !!   dia_dct_wri  :  Write tranports results in ascii files 
     24   !!   interp       :  Compute temperature/salinity/density at U-point or V-point 
     25   !!    
     26   !!---------------------------------------------------------------------- 
     27   USE oce             ! ocean dynamics and tracers 
     28   USE dom_oce         ! ocean space and time domain 
     29   USE phycst          ! physical constants 
     30   USE in_out_manager  ! I/O manager 
     31   USE daymod          ! calendar 
     32   USE dianam          ! build name of file 
     33   USE lib_mpp         ! distributed memory computing library 
    4034#if defined key_lim2 
    41   USE ice_2 
     35   USE ice_2 
    4236#endif 
    4337#if defined key_lim3 
    44   USE ice 
     38   USE ice 
    4539#endif 
    46   USE domvvl 
    47   USE timing          ! preformance summary 
    48   USE wrk_nemo        ! working arrays 
    49  
    50   IMPLICIT NONE 
    51   PRIVATE 
    52  
    53   !! * Routine accessibility 
    54   PUBLIC   dia_dct      ! routine called by step.F90 
    55   PUBLIC   dia_dct_init ! routine called by opa.F90 
    56   PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90  
    57   PRIVATE  readsec 
    58   PRIVATE  removepoints 
    59   PRIVATE  transport 
    60   PRIVATE  dia_dct_wri 
    61  
    62 #include "domzgr_substitute.h90" 
    63  
    64   !! * Shared module variables 
    65   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag 
    66  
    67   !! * Module variables 
    68   INTEGER :: nn_dct        ! Frequency of computation 
    69   INTEGER :: nn_dctwri     ! Frequency of output 
    70   INTEGER :: nn_secdebug   ! Number of the section to debug 
     40   USE domvvl 
     41   USE timing          ! preformance summary 
     42   USE wrk_nemo        ! working arrays 
     43 
     44   IMPLICIT NONE 
     45   PRIVATE 
     46 
     47   PUBLIC   dia_dct      ! routine called by step.F90 
     48   PUBLIC   dia_dct_init ! routine called by opa.F90 
     49   PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90  
     50   PRIVATE  readsec 
     51   PRIVATE  removepoints 
     52   PRIVATE  transport 
     53   PRIVATE  dia_dct_wri 
     54 
     55   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag 
     56 
     57   INTEGER :: nn_dct        ! Frequency of computation 
     58   INTEGER :: nn_dctwri     ! Frequency of output 
     59   INTEGER :: nn_secdebug   ! Number of the section to debug 
    7160    
    72   INTEGER, PARAMETER :: nb_class_max  = 10 
    73   INTEGER, PARAMETER :: nb_sec_max    = 150 
    74   INTEGER, PARAMETER :: nb_point_max  = 2000 
    75   INTEGER, PARAMETER :: nb_type_class = 10 
    76   INTEGER, PARAMETER :: nb_3d_vars    = 3  
    77   INTEGER, PARAMETER :: nb_2d_vars    = 2  
    78   INTEGER            :: nb_sec  
    79  
    80   TYPE POINT_SECTION 
    81      INTEGER :: I,J 
    82   END TYPE POINT_SECTION 
    83  
    84   TYPE COORD_SECTION 
    85      REAL(wp) :: lon,lat 
    86   END TYPE COORD_SECTION 
    87  
    88   TYPE SECTION 
    89      CHARACTER(len=60)                            :: name              ! name of the sec 
    90      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and 
     61   INTEGER, PARAMETER :: nb_class_max  = 10 
     62   INTEGER, PARAMETER :: nb_sec_max    = 150 
     63   INTEGER, PARAMETER :: nb_point_max  = 2000 
     64   INTEGER, PARAMETER :: nb_type_class = 10 
     65   INTEGER, PARAMETER :: nb_3d_vars    = 3  
     66   INTEGER, PARAMETER :: nb_2d_vars    = 2  
     67   INTEGER            :: nb_sec  
     68 
     69   TYPE POINT_SECTION 
     70      INTEGER :: I,J 
     71   END TYPE POINT_SECTION 
     72 
     73   TYPE COORD_SECTION 
     74      REAL(wp) :: lon,lat 
     75   END TYPE COORD_SECTION 
     76 
     77   TYPE SECTION 
     78      CHARACTER(len=60)                            :: name              ! name of the sec 
     79      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and 
    9180                                                                       ! heat transports 
    92      LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation 
    93      LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line 
    94      TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec 
    95      INTEGER                                      :: nb_class          ! number of boundaries for density classes 
    96      INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section 
    97      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class 
    98      REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want) 
    99                                                      zsigp           ,&! potential density classes    (99 if you don't want) 
    100                                                      zsal            ,&! salinity classes   (99 if you don't want) 
    101                                                      ztem            ,&! temperature classes(99 if you don't want) 
    102                                                      zlay              ! level classes      (99 if you don't want) 
    103      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
    104      REAL(wp)                                         :: slopeSection  ! slope of the section 
    105      INTEGER                                          :: nb_point      ! number of points in the section 
    106      TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections 
    107   END TYPE SECTION 
    108  
    109   TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections 
    110   
    111   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d  
    112   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
    113  
     81      LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation 
     82      LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line 
     83      TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec 
     84      INTEGER                                      :: nb_class          ! number of boundaries for density classes 
     85      INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section 
     86      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class 
     87      REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want) 
     88                                                      zsigp           ,&! potential density classes    (99 if you don't want) 
     89                                                      zsal            ,&! salinity classes   (99 if you don't want) 
     90                                                      ztem            ,&! temperature classes(99 if you don't want) 
     91                                                      zlay              ! level classes      (99 if you don't want) 
     92      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
     93      REAL(wp)                                         :: slopeSection  ! slope of the section 
     94      INTEGER                                          :: nb_point      ! number of points in the section 
     95      TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections 
     96   END TYPE SECTION 
     97 
     98   TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections 
     99  
     100   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d  
     101   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
     102 
     103   !!---------------------------------------------------------------------- 
     104   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    114105   !! $Id$ 
     106   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     107   !!---------------------------------------------------------------------- 
    115108CONTAINS 
    116  
    117109  
    118110  INTEGER FUNCTION diadct_alloc()  
     
    130122  
    131123  END FUNCTION diadct_alloc  
     124 
    132125 
    133126  SUBROUTINE dia_dct_init 
     
    208201     !!               Reinitialise all relevant arrays to zero  
    209202     !!--------------------------------------------------------------------- 
    210      !! * Arguments 
    211      INTEGER,INTENT(IN)        ::kt 
    212  
    213      !! * Local variables 
     203     INTEGER,INTENT(in)        ::kt 
     204     ! 
    214205     INTEGER             :: jsec,            &! loop on sections 
    215206                            itotal            ! nb_sec_max*nb_type_class*nb_class_max 
     
    220211     REAL(wp), POINTER, DIMENSION(:)    :: zwork !   "   
    221212     REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   " 
    222  
    223213     !!---------------------------------------------------------------------     
     214     ! 
    224215     IF( nn_timing == 1 )   CALL timing_start('dia_dct') 
    225216 
     
    619610                            zumid_ice, zvmid_ice,                &!U/V ice velocity  
    620611                            zTnorm                                !transport of velocity through one cell's sides  
    621      REAL(wp)            :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point 
     612     REAL(wp)            :: ztn, zsn, zrhoi, zrhop, zsshn, zdep !temperature/salinity/potential density/ssh/depth at u/v point 
    622613 
    623614     TYPE(POINT_SECTION) :: k 
    624      !!-------------------------------------------------------- 
    625  
    626      IF( ld_debug )WRITE(numout,*)'      Compute transport' 
    627  
    628      !---------------------------! 
    629      !  COMPUTE TRANSPORT        ! 
    630      !---------------------------! 
    631      IF(sec%nb_point .NE. 0)THEN    
    632  
    633         !---------------------------------------------------------------------------------------------------- 
    634         !Compute sign for velocities: 
    635         ! 
    636         !convention: 
    637         !   non horizontal section: direction + is toward left hand of section 
    638         !       horizontal section: direction + is toward north of section 
    639         ! 
    640         ! 
    641         !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0 
    642         !       ----------------      -----------------     ---------------             -------------- 
    643         ! 
    644         !   isgnv=1         direction +       
    645         !  ______         _____             ______                                                    
    646         !        |           //|            |                  |                         direction +    
    647         !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\ 
    648         !        |_______  //         ______|    \\            | ---\                        | 
    649         !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________ 
    650         !               |             |          __\\|         |                     
    651         !               |             |     direction +        |                      isgnv=1                                  
    652         !                                                       
    653         !---------------------------------------------------------------------------------------------------- 
    654         isgnu = 1 
    655         IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1  
    656         ELSE                                ; isgnv =  1 
    657         ENDIF 
    658         IF( sec%slopeSection .GE. 9999. )     isgnv =  1 
    659  
    660         IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv 
    661  
    662         !--------------------------------------! 
    663         ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
    664         !--------------------------------------! 
    665         DO jseg=1,MAX(sec%nb_point-1,0) 
     615      !!-------------------------------------------------------- 
     616      ! 
     617      IF( ld_debug )WRITE(numout,*)'      Compute transport' 
     618 
     619      !---------------------------! 
     620      !  COMPUTE TRANSPORT        ! 
     621      !---------------------------! 
     622      IF(sec%nb_point .NE. 0)THEN    
     623 
     624         !---------------------------------------------------------------------------------------------------- 
     625         !Compute sign for velocities: 
     626         ! 
     627         !convention: 
     628         !   non horizontal section: direction + is toward left hand of section 
     629         !       horizontal section: direction + is toward north of section 
     630         ! 
     631         ! 
     632         !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0 
     633         !       ----------------      -----------------     ---------------             -------------- 
     634         ! 
     635         !   isgnv=1         direction +       
     636         !  ______         _____             ______                                                    
     637         !        |           //|            |                  |                         direction +    
     638         !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\ 
     639         !        |_______  //         ______|    \\            | ---\                        | 
     640         !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________ 
     641         !               |             |          __\\|         |                     
     642         !               |             |     direction +        |                      isgnv=1                                  
     643         !                                                       
     644         !---------------------------------------------------------------------------------------------------- 
     645         isgnu = 1 
     646         IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1  
     647         ELSE                                ; isgnv =  1 
     648         ENDIF 
     649         IF( sec%slopeSection .GE. 9999. )     isgnv =  1 
     650 
     651         IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv 
     652 
     653         !--------------------------------------! 
     654         ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     655         !--------------------------------------! 
     656         DO jseg=1,MAX(sec%nb_point-1,0) 
    666657               
    667            !------------------------------------------------------------------------------------------- 
    668            ! Select the appropriate coordinate for computing the velocity of the segment 
    669            ! 
    670            !                      CASE(0)                                    Case (2) 
    671            !                      -------                                    -------- 
    672            !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
    673            !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
    674            !                                                                            | 
    675            !                                                                            | 
    676            !                                                                            | 
    677            !                      Case (3)                                            U(i,j) 
    678            !                      --------                                              | 
    679            !                                                                            | 
    680            !  listPoint(jseg+1) F(i,j+1)                                                | 
    681            !                        |                                                   | 
    682            !                        |                                                   | 
    683            !                        |                                 listPoint(jseg+1) F(i,j-1) 
    684            !                        |                                             
    685            !                        |                                             
    686            !                     U(i,j+1)                                             
    687            !                        |                                       Case(1)      
    688            !                        |                                       ------       
    689            !                        |                                             
    690            !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
    691            !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
    692            ! listPoint(jseg)     F(i,j) 
    693            !  
    694            !------------------------------------------------------------------------------------------- 
    695  
    696            SELECT CASE( sec%direction(jseg) ) 
    697            CASE(0)  ;   k = sec%listPoint(jseg) 
    698            CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
    699            CASE(2)  ;   k = sec%listPoint(jseg) 
    700            CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
    701            END SELECT 
    702  
    703            !---------------------------|  
    704            !     LOOP ON THE LEVEL     |  
    705            !---------------------------|  
    706            !Sum of the transport on the vertical   
    707            DO jk=1,mbathy(k%I,k%J)  
    708   
    709               ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point  
    710               SELECT CASE( sec%direction(jseg) )  
    711               CASE(0,1)  
    712                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
    713                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    714                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
    715                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
    716                  zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)  
    717               CASE(2,3)  
    718                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
    719                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    720                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
    721                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
    722                  zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    723               END SELECT  
    724   
    725               zfsdep= fsdept(k%I,k%J,jk)  
     658            !------------------------------------------------------------------------------------------- 
     659            ! Select the appropriate coordinate for computing the velocity of the segment 
     660            ! 
     661            !                      CASE(0)                                    Case (2) 
     662            !                      -------                                    -------- 
     663            !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     664            !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     665            !                                                                            | 
     666            !                                                                            | 
     667            !                                                                            | 
     668            !                      Case (3)                                            U(i,j) 
     669            !                      --------                                              | 
     670            !                                                                            | 
     671            !  listPoint(jseg+1) F(i,j+1)                                                | 
     672            !                        |                                                   | 
     673            !                        |                                                   | 
     674            !                        |                                 listPoint(jseg+1) F(i,j-1) 
     675            !                        |                                             
     676            !                        |                                             
     677            !                     U(i,j+1)                                             
     678            !                        |                                       Case(1)      
     679            !                        |                                       ------       
     680            !                        |                                             
     681            !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     682            !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     683            ! listPoint(jseg)     F(i,j) 
     684            !  
     685            !------------------------------------------------------------------------------------------- 
     686 
     687            SELECT CASE( sec%direction(jseg) ) 
     688            CASE(0)   ;    k = sec%listPoint(jseg) 
     689            CASE(1)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     690            CASE(2)   ;    k = sec%listPoint(jseg) 
     691            CASE(3)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     692            END SELECT 
     693 
     694            !---------------------------|  
     695            !     LOOP ON THE LEVEL     |  
     696            !---------------------------|  
     697            DO jk = 1, mbathy(k%I,k%J)            !Sum of the transport on the vertical 
     698            !           ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point  
     699            SELECT CASE( sec%direction(jseg) ) 
     700               CASE(0,1)  
     701                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
     702                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
     703                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
     704                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
     705                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)  
     706               CASE(2,3)  
     707                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
     708                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
     709                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
     710                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
     711                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
     712               END SELECT  
     713               ! 
     714               zdep= gdept_n(k%I,k%J,jk)  
    726715   
    727               !compute velocity with the correct direction  
    728               SELECT CASE( sec%direction(jseg) )  
    729               CASE(0,1)    
    730                  zumid=0.  
    731                  zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)  
    732               CASE(2,3)  
    733                  zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)  
    734                  zvmid=0.  
    735               END SELECT  
    736   
    737               !zTnorm=transport through one cell;  
    738               !velocity* cell's length * cell's thickness  
    739               zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &  
    740                      zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)  
    741  
    742 #if ! defined key_vvl 
    743               !add transport due to free surface  
    744               IF( jk==1 )THEN  
    745                  zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &  
    746                                    zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)  
    747               ENDIF  
    748 #endif 
     716               SELECT CASE( sec%direction(jseg) )                !compute velocity with the correct direction  
     717               CASE(0,1)   
     718                  zumid=0._wp 
     719                  zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)  
     720               CASE(2,3)  
     721                  zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)  
     722                  zvmid=0._wp 
     723               END SELECT  
     724  
     725               !zTnorm=transport through one cell;  
     726               !velocity* cell's length * cell's thickness  
     727               zTnorm = zumid*e2u(k%I,k%J) * e3u_n(k%I,k%J,jk)     &  
     728                  &   + zvmid*e1v(k%I,k%J) * e3v_n(k%I,k%J,jk)  
     729 
     730!!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!! 
     731               IF( ln_linssh ) THEN              !add transport due to free surface  
     732                  IF( jk==1 ) THEN  
     733                     zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk)   &  
     734                        &            + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)  
     735                  ENDIF  
     736               ENDIF 
     737!!gm end 
    749738              !COMPUTE TRANSPORT   
    750739  
    751740              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm  
    752741   
    753               IF ( sec%llstrpond ) THEN  
     742              IF( sec%llstrpond ) THEN  
    754743                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp 
    755744                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001 
    756745              ENDIF 
    757746    
    758            ENDDO !end of loop on the level 
     747           END DO !end of loop on the level 
    759748 
    760749#if defined key_lim2 || defined key_lim3 
     
    797786                 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   & 
    798787                                   a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) 
    799               ENDDO 
     788              END DO 
    800789#endif 
    801790    
     
    803792#endif 
    804793  
    805         ENDDO !end of loop on the segment 
     794        END DO !end of loop on the segment 
    806795 
    807796     ENDIF !end of sec%nb_point =0 case 
    808797     ! 
    809798  END SUBROUTINE transport 
    810    
     799 
     800 
    811801  SUBROUTINE dia_dct_sum(sec,jsec)  
    812802     !!-------------------------------------------------------------  
     
    828818     !!  
    829819     !!-------------------------------------------------------------  
    830      !! * arguments  
    831820     TYPE(SECTION),INTENT(INOUT) :: sec  
    832821     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section  
     
    834823     TYPE(POINT_SECTION) :: k  
    835824     INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes   
    836      REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point  
     825     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point  
    837826     !!-------------------------------------------------------------  
    838827  
     
    903892              END SELECT  
    904893  
    905               zfsdep= fsdept(k%I,k%J,jk)  
     894              zdep= gdept_n(k%I,k%J,jk)  
    906895   
    907896              !-------------------------------  
     
    932921                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     &  
    933922  
    934                     ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &  
    935                     (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &  
     923                    ((( zdep .GE. sec%zlay(jclass)) .AND.                &  
     924                    (   zdep .LE. sec%zlay(jclass+1))) .OR.              &  
    936925                    ( sec%zlay(jclass) .EQ. 99. ))                         &  
    937926                                                                   ))   THEN  
     
    991980#endif  
    992981   
    993         ENDDO !end of loop on the segment  
     982        END DO !end of loop on the segment  
    994983  
    995984     ELSE  !if sec%nb_point =0  
     
    1000989  
    1001990  END SUBROUTINE dia_dct_sum  
    1002    
     991 
     992 
    1003993  SUBROUTINE dia_dct_wri(kt,ksec,sec) 
    1004994     !!------------------------------------------------------------- 
     
    11381128                              sec%transport(9,1),sec%transport(10,1), & 
    11391129                              sec%transport(9,1)+sec%transport(10,1)  
    1140      ENDIF 
     1130      ENDIF 
    11411131                                               
    1142 118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4) 
    1143 119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 
    1144  
    1145      CALL wrk_dealloc(nb_type_class , zsumclasses )   
    1146   END SUBROUTINE dia_dct_wri 
    1147  
    1148   FUNCTION interp(ki, kj, kk, cd_point, ptab) 
     1132118   FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4) 
     1133119   FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 
     1134 
     1135      CALL wrk_dealloc(nb_type_class , zsumclasses )   
     1136      ! 
     1137   END SUBROUTINE dia_dct_wri 
     1138 
     1139 
     1140   FUNCTION interp(ki, kj, kk, cd_point, ptab) 
    11491141  !!---------------------------------------------------------------------- 
    11501142  !! 
     
    12141206  !*local declations 
    12151207  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer 
    1216   REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real 
     1208  REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu            ! local real 
    12171209  REAL(wp):: zet1, zet2                                        ! weight for interpolation  
    12181210  REAL(wp):: zdep1,zdep2                                       ! differences of depth 
     
    12411233  IF( ln_sco )THEN   ! s-coordinate case 
    12421234 
    1243      zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2  
    1244      zdep1 = fsdept(ii1,ij1,kk) - zdepu 
    1245      zdep2 = fsdept(ii2,ij2,kk) - zdepu 
     1235     zdepu = ( gdept_n(ii1,ij1,kk) +  gdept_n(ii2,ij2,kk) ) * 0.5_wp  
     1236     zdep1 = gdept_n(ii1,ij1,kk) - zdepu 
     1237     zdep2 = gdept_n(ii2,ij2,kk) - zdepu 
    12461238 
    12471239     ! weights 
     
    12551247  ELSE       ! full step or partial step case  
    12561248 
    1257 #if defined key_vvl 
    1258  
    1259      ze3t  = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk)  
    1260      zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk) 
    1261      zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk) 
    1262  
    1263 #else 
    1264  
    1265      ze3t  = fse3t(ii2,ij2,kk)   - fse3t(ii1,ij1,kk)  
    1266      zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk) 
    1267      zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk) 
    1268  
    1269 #endif 
     1249     ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk)  
     1250     zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk) 
     1251     zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk) 
    12701252 
    12711253     IF(kk .NE. 1)THEN 
     
    12881270 
    12891271  ENDIF 
    1290  
    1291  
    1292   END FUNCTION interp 
     1272      ! 
     1273   END FUNCTION interp 
    12931274 
    12941275#else 
     
    13111292#endif 
    13121293 
     1294   !!====================================================================== 
    13131295END MODULE diadct 
Note: See TracChangeset for help on using the changeset viewer.