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 11134 for branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 – NEMO

Ignore:
Timestamp:
2019-06-18T17:48:39+02:00 (5 years ago)
Author:
jcastill
Message:

Full set of changes as in the original branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    r11132 r11134  
    5959  PRIVATE  transport 
    6060  PRIVATE  dia_dct_wri 
     61  PRIVATE  dia_dct_wri_NOOS 
    6162 
    6263#include "domzgr_substitute.h90" 
     
    6667 
    6768  !! * 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 
     69  INTEGER :: nn_dct      = 1     ! Frequency of computation 
     70  INTEGER :: nn_dctwri   = 1     ! Frequency of output 
     71  INTEGER :: nn_secdebug = 0     ! Number of the section to debug 
     72  INTEGER :: nn_dct_h    = 1     ! Frequency of computation for NOOS hourly files 
     73  INTEGER :: nn_dctwri_h = 1     ! Frequency of output for NOOS hourly files 
    7174    
    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  
     75  INTEGER, PARAMETER :: nb_class_max  = 12   ! maximum number of classes, i.e. depth levels or density classes 
     76  INTEGER, PARAMETER :: nb_sec_max    = 30   ! maximum number of sections 
     77  INTEGER, PARAMETER :: nb_point_max  = 300  ! maximum number of points in a single section 
     78  INTEGER, PARAMETER :: nb_type_class       = 14   ! types of calculations, i.e. pos transport, neg transport, heat transport, salt transport 
     79  INTEGER, PARAMETER :: nb_3d_vars    = 5 
     80  INTEGER, PARAMETER :: nb_2d_vars    = 2 
    7881  INTEGER            :: nb_sec  
    7982 
     
    102105                                                     zlay              ! level classes      (99 if you don't want) 
    103106     REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
     107     REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport_h   ! transport output 
    104108     REAL(wp)                                         :: slopeSection  ! slope of the section 
    105109     INTEGER                                          :: nb_point      ! number of points in the section 
     
    111115  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d  
    112116  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
     117  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d_h 
     118  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d_h 
     119  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  z_hr_output 
    113120 
    114121   !! $Id$ 
     
    120127     !!                   ***  FUNCTION diadct_alloc  ***  
    121128     !!----------------------------------------------------------------------  
    122      INTEGER :: ierr(2)  
     129     INTEGER :: ierr(5)  
    123130     !!----------------------------------------------------------------------  
    124131 
    125      ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) )  
    126      ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) )  
     132     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk)  , STAT=ierr(1) )  
     133     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)      , STAT=ierr(2) )  
     134     ALLOCATE(transports_3d_h(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(3) ) 
     135     ALLOCATE(transports_2d_h(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(4) ) 
     136     ALLOCATE(z_hr_output(nb_sec_max,168,nb_class_max)               , STAT=ierr(5) ) ! 168 = 24 * 7days 
    127137 
    128138     diadct_alloc = MAXVAL( ierr )  
     
    153163     IF(lwm) WRITE ( numond, namdct ) 
    154164 
     165     IF( ln_NOOS ) THEN 
     166        nn_dct=3600./rdt         ! hard coded for NOOS transects, to give 25 hour means  
     167        nn_dctwri=86400./rdt 
     168        nn_dct_h=1       ! hard coded for NOOS transects, to give hourly data 
     169        nn_dctwri_h=3600./rdt 
     170     ENDIF 
     171 
    155172     IF( lwp ) THEN 
    156173        WRITE(numout,*) " " 
    157174        WRITE(numout,*) "diadct_init: compute transports through sections " 
    158175        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
    159         WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
    160         WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     176        IF( ln_NOOS ) THEN 
     177           WRITE(numout,*) "       Frequency of computation hard coded to be every hour: nn_dct    = ",nn_dct 
     178           WRITE(numout,*) "       Frequency of write hard coded to average 25 instantaneous hour values: nn_dctwri = ",nn_dctwri 
     179           WRITE(numout,*) "       Frequency of hourly computation hard coded to be every timestep: nn_dct_h  = ",nn_dct_h 
     180           WRITE(numout,*) "       Frequency of hourly write hard coded to every hour: nn_dctwri_h = ",nn_dctwri_h 
     181        ELSE 
     182           WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
     183           WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     184        ENDIF 
    161185 
    162186        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN 
     
    177201     !open output file 
    178202     IF( lwm ) THEN 
    179         CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    180         CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    181         CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     203        IF( ln_NOOS ) THEN 
     204           CALL ctl_opn( numdct_NOOS  ,'NOOS_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     205           CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     206        ELSE 
     207           CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     208           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     209           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     210        ENDIF 
    182211     ENDIF 
    183212 
    184213     ! Initialise arrays to zero  
    185      transports_3d(:,:,:,:)=0.0  
    186      transports_2d(:,:,:)  =0.0  
     214     transports_3d(:,:,:,:)  =0.0  
     215     transports_2d(:,:,:)    =0.0  
     216     transports_3d_h(:,:,:,:)=0._wp 
     217     transports_2d_h(:,:,:)  =0._wp 
     218     z_hr_output(:,:,:)      =0._wp 
    187219 
    188220     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init') 
     
    229261        CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum  ) 
    230262     ENDIF     
    231   
     263     lldebug=.TRUE. 
    232264     ! Initialise arrays 
    233265     zwork(:) = 0.0  
     
    243275  
    244276     ! Compute transport and write only at nn_dctwri 
    245      IF( MOD(kt,nn_dct)==0 ) THEN  
     277     IF( MOD(kt,nn_dct)==0 .or. &                ! compute transport every nn_dct time steps 
     278         (ln_NOOS .and. kt==nn_it000 ) )  THEN   ! also include first time step when calculating NOOS 25 hour averages 
    246279 
    247280        DO jsec=1,nb_sec 
    248281 
    249282           !debug this section computing ? 
    250            lldebug=.FALSE. 
    251283           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE.  
    252284 
     
    271303           !Sum on all procs  
    272304           IF( lk_mpp )THEN 
     305              zsum(:,:,:)=0.0_wp 
    273306              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max  
    274307              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/) 
     
    283316           DO jsec=1,nb_sec 
    284317 
    285               IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     318              IF( lwm .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     319              IF( lwm .and.       ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting 
    286320             
    287321              !nullify transports values after writing 
    288322              transports_3d(:,jsec,:,:)=0. 
    289323              transports_2d(:,jsec,:  )=0. 
    290               secs(jsec)%transport(:,:)=0.   
     324              secs(jsec)%transport(:,:)=0.  
     325              IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values) 
    291326 
    292327           ENDDO 
     
    295330 
    296331     ENDIF 
     332 
     333     IF ( MOD(kt,nn_dct_h)==0 ) THEN            ! compute transport every nn_dct_h time steps 
     334 
     335        DO jsec=1,nb_sec 
     336 
     337           !lldebug=.FALSE. 
     338           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE.  
     339 
     340           !Compute transport through section   
     341           CALL transport_h(secs(jsec),lldebug,jsec)  
     342 
     343        ENDDO 
     344              
     345        IF( MOD(kt,nn_dctwri_h)==0 )THEN 
     346 
     347           IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)"      diadct: average and write hourly files at kt = ",kt          
     348   
     349           !! divide arrays by nn_dctwri/nn_dct to obtain average 
     350           transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h) 
     351           transports_2d_h(:,:,:)  =transports_2d_h(:,:,:)  /(nn_dctwri_h/nn_dct_h) 
     352 
     353           ! Sum over each class 
     354           DO jsec=1,nb_sec 
     355              CALL dia_dct_sum_h(secs(jsec),jsec) 
     356           ENDDO 
     357  
     358           !Sum on all procs  
     359          IF( lk_mpp )THEN 
     360              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max  
     361              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/) 
     362              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO 
     363              zwork(:)= RESHAPE(zsum(:,:,:), ish ) 
     364              CALL mpp_sum(zwork, ish(1)) 
     365              zsum(:,:,:)= RESHAPE(zwork,ish2) 
     366              DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO 
     367           ENDIF 
     368 
     369           !Write the transport 
     370           DO jsec=1,nb_sec 
     371 
     372              IF( lwp .and.       ln_NOOS )CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec))   ! use NOOS specific formatting 
     373             
     374              !nullify transports values after writing 
     375              transports_3d_h(:,jsec,:,:)=0.0 
     376              transports_2d_h(:,jsec,:)=0.0 
     377              secs(jsec)%transport_h(:,:)=0.   
     378              IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values) 
     379 
     380           ENDDO 
     381 
     382        ENDIF  
     383 
     384     ENDIF     
    297385 
    298386     IF( lk_mpp )THEN 
     
    356444        secs(jsec)%zlay         = 99._wp          
    357445        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0 
     446        secs(jsec)%transport_h  =  0._wp   ; secs(jsec)%nb_point       = 0 
    358447 
    359448        !read section's number / name / computing choices / classes / slopeSection / points number 
    360449        !----------------------------------------------------------------------------------------- 
    361450        READ(numdct_in,iostat=iost)isec 
    362         IF (iost .NE. 0 )EXIT !end of file  
     451        IF (iost .NE. 0 )EXIT !end of file 
    363452        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec 
    364         IF( jsec .NE. isec )  CALL ctl_stop( cltmp ) 
     453        IF( jsec .NE. isec )CALL ctl_stop( cltmp ) 
    365454 
    366455        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec  
     
    380469        READ(numdct_in)iptglo 
    381470 
     471        IF ( ln_NOOS ) THEN 
     472           WRITE(numout,*) 'Section name = ',TRIM(secs(jsec)%name) 
     473           WRITE(numout,*) 'coordSec = ',secs(jsec)%coordSec 
     474           WRITE(numout,*) 'iptglo = ',iptglo 
     475        ENDIF 
     476 
    382477        !debug 
    383478        !----- 
     
    405500           !read points'coordinates and directions  
    406501           !-------------------------------------- 
     502           IF ( ln_NOOS ) THEN 
     503              WRITE(numout,*) 'Coords and direction read' 
     504           ENDIF 
     505 
    407506           coordtemp(:) = POINT_SECTION(0,0) !list of points read 
    408507           directemp(:) = 0                  !value of directions of each points 
     
    422521              ENDDO                   
    423522           ENDIF 
    424   
     523 
    425524           !Now each proc selects only points that are in its domain: 
    426525           !-------------------------------------------------------- 
     
    624723     !!-------------------------------------------------------- 
    625724 
    626      IF( ld_debug )WRITE(numout,*)'      Compute transport' 
     725     !!NIALL IF( ld_debug )WRITE(numout,*)'      Compute transport' 
    627726 
    628727     !---------------------------! 
     
    710809              SELECT CASE( sec%direction(jseg) )  
    711810              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)  
     811                 ztn   = interp(k%I,k%J,jk,'V',0 )  
     812                 zsn   = interp(k%I,k%J,jk,'V',1 )  
     813                 zrhop = interp(k%I,k%J,jk,'V',2 )  
     814                 zrhoi = interp(k%I,k%J,jk,'V',3 )  
    716815                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)  
    717816              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)  
     817                 ztn   = interp(k%I,k%J,jk,'U',0 )  
     818                 zsn   = interp(k%I,k%J,jk,'U',1 )  
     819                 zrhop = interp(k%I,k%J,jk,'U',2 )  
     820                 zrhoi = interp(k%I,k%J,jk,'U',3 )  
    722821                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    723822              END SELECT  
     
    754853                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp 
    755854                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001 
     855                 transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 
     856                 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp 
    756857              ENDIF 
    757858    
     
    809910  END SUBROUTINE transport 
    810911   
     912  SUBROUTINE transport_h(sec,ld_debug,jsec) 
     913     !!------------------------------------------------------------------------------------------- 
     914     !!                     ***  ROUTINE hourly_transport  *** 
     915     !! 
     916     !!              Exactly as routine transport but for data summed at 
     917     !!              each time step and averaged each hour 
     918     !! 
     919     !!  Purpose ::  Compute the transport for each point in a section 
     920     !! 
     921     !!  Method  ::  Loop over each segment, and each vertical level and add the transport 
     922     !!              Be aware :           
     923     !!              One section is a sum of segments 
     924     !!              One segment is defined by 2 consecutive points in sec%listPoint 
     925     !!              All points of sec%listPoint are positioned on the F-point of the cell 
     926     !!  
     927     !!              There are two loops:                  
     928     !!              loop on the segment between 2 nodes 
     929     !!              loop on the level jk 
     930     !! 
     931     !! ** Output: Arrays containing the volume,density,salinity,temperature etc 
     932     !!            transports for each point in a section, summed over each nn_dct. 
     933     !! 
     934     !!------------------------------------------------------------------------------------------- 
     935     !! * Arguments 
     936     TYPE(SECTION),INTENT(INOUT) :: sec 
     937     LOGICAL      ,INTENT(IN)    :: ld_debug 
     938     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     939     
     940     !! * Local variables 
     941     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes  
     942                            isgnu  , isgnv      ! 
     943     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment 
     944                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity 
     945                zTnorm                          !transport of velocity through one cell's sides 
     946     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     947 
     948     TYPE(POINT_SECTION) :: k 
     949     !!-------------------------------------------------------- 
     950 
     951     !!NIALL IF( ld_debug )WRITE(numout,*)'      Compute transport' 
     952 
     953     !---------------------------! 
     954     !  COMPUTE TRANSPORT        ! 
     955     !---------------------------! 
     956     IF(sec%nb_point .NE. 0)THEN    
     957 
     958        !---------------------------------------------------------------------------------------------------- 
     959        !Compute sign for velocities: 
     960        ! 
     961        !convention: 
     962        !   non horizontal section: direction + is toward left hand of section 
     963        !       horizontal section: direction + is toward north of section 
     964        ! 
     965        ! 
     966        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0 
     967        !       ----------------      -----------------     ---------------             -------------- 
     968        ! 
     969        !   isgnv=1         direction +       
     970        !  ______         _____             ______                                                    
     971        !        |           //|            |                  |                         direction +    
     972        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\ 
     973        !        |_______  //         ______|    \\            | ---\                        | 
     974        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________ 
     975        !               |             |          __\\|         |                     
     976        !               |             |     direction +        |                      isgnv=1                                  
     977        !                                                       
     978        !---------------------------------------------------------------------------------------------------- 
     979        isgnu = 1 
     980        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1  
     981        ELSE                                ; isgnv =  1 
     982        ENDIF 
     983 
     984        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv 
     985 
     986        !--------------------------------------! 
     987        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     988        !--------------------------------------! 
     989        DO jseg=1,MAX(sec%nb_point-1,0) 
     990               
     991           !------------------------------------------------------------------------------------------- 
     992           ! Select the appropriate coordinate for computing the velocity of the segment 
     993           ! 
     994           !                      CASE(0)                                    Case (2) 
     995           !                      -------                                    -------- 
     996           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     997           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     998           !                                                                            | 
     999           !                                                                            | 
     1000           !                                                                            | 
     1001           !                      Case (3)                                            U(i,j) 
     1002           !                      --------                                              | 
     1003           !                                                                            | 
     1004           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1005           !                        |                                                   | 
     1006           !                        |                                                   | 
     1007           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1008           !                        |                                             
     1009           !                        |                                             
     1010           !                     U(i,j+1)                                             
     1011           !                        |                                       Case(1)      
     1012           !                        |                                       ------       
     1013           !                        |                                             
     1014           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1015           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1016           ! listPoint(jseg)     F(i,j) 
     1017           !  
     1018           !------------------------------------------------------------------------------------------- 
     1019 
     1020           SELECT CASE( sec%direction(jseg) ) 
     1021           CASE(0)  ;   k = sec%listPoint(jseg) 
     1022           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1023           CASE(2)  ;   k = sec%listPoint(jseg) 
     1024           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1025           END SELECT 
     1026 
     1027           !---------------------------| 
     1028           !     LOOP ON THE LEVEL     | 
     1029           !---------------------------| 
     1030           !Sum of the transport on the vertical  
     1031           DO jk=1,mbathy(k%I,k%J) 
     1032 
     1033              ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 
     1034              SELECT CASE( sec%direction(jseg) ) 
     1035              CASE(0,1) 
     1036                 ztn   = interp(k%I,k%J,jk,'V',0 ) 
     1037                 zsn   = interp(k%I,k%J,jk,'V',1) 
     1038                 zrhop = interp(k%I,k%J,jk,'V',2) 
     1039                 zrhoi = interp(k%I,k%J,jk,'V',3) 
     1040                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1041              CASE(2,3) 
     1042                 ztn   = interp(k%I,k%J,jk,'U',0) 
     1043                 zsn   = interp(k%I,k%J,jk,'U',1) 
     1044                 zrhop = interp(k%I,k%J,jk,'U',2) 
     1045                 zrhoi = interp(k%I,k%J,jk,'U',3) 
     1046                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1047              END SELECT 
     1048 
     1049              zfsdep= fsdept(k%I,k%J,jk) 
     1050  
     1051              !compute velocity with the correct direction 
     1052              SELECT CASE( sec%direction(jseg) ) 
     1053              CASE(0,1)   
     1054                 zumid=0. 
     1055                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
     1056              CASE(2,3) 
     1057                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
     1058                 zvmid=0. 
     1059              END SELECT 
     1060 
     1061              !zTnorm=transport through one cell; 
     1062              !velocity* cell's length * cell's thickness 
     1063              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
     1064                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
     1065 
     1066#if ! defined key_vvl 
     1067              !add transport due to free surface 
     1068              IF( jk==1 )THEN 
     1069                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
     1070                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
     1071              ENDIF 
     1072#endif 
     1073              !COMPUTE TRANSPORT  
     1074 
     1075              transports_3d_h(1,jsec,jseg,jk) = transports_3d_h(1,jsec,jseg,jk) + zTnorm 
     1076  
     1077              IF ( sec%llstrpond ) THEN 
     1078                 transports_3d_h(2,jsec,jseg,jk) = transports_3d_h(2,jsec,jseg,jk)  + zTnorm * zrhoi 
     1079                 transports_3d_h(3,jsec,jseg,jk) = transports_3d_h(3,jsec,jseg,jk)  + zTnorm * zrhop 
     1080                 transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 
     1081                 transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp 
     1082              ENDIF 
     1083 
     1084           ENDDO !end of loop on the level 
     1085 
     1086#if defined key_lim2 || defined key_lim3 
     1087 
     1088           !ICE CASE     
     1089           !------------ 
     1090           IF( sec%ll_ice_section )THEN 
     1091              SELECT CASE (sec%direction(jseg)) 
     1092              CASE(0) 
     1093                 zumid_ice = 0 
     1094                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 
     1095              CASE(1) 
     1096                 zumid_ice = 0 
     1097                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 
     1098              CASE(2) 
     1099                 zvmid_ice = 0 
     1100                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 
     1101              CASE(3) 
     1102                 zvmid_ice = 0 
     1103                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 
     1104              END SELECT 
     1105    
     1106              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 
     1107    
     1108              transports_2d_h(1,jsec,jseg) = transports_2d_h(1,jsec,jseg) + (zTnorm)*   & 
     1109                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1110                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
     1111                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 
     1112                                   +zice_vol_pos 
     1113              transports_2d_h(2,jsec,jseg) = transports_2d_h(2,jsec,jseg) + (zTnorm)*   & 
     1114                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1115                                   +zice_surf_pos 
     1116    
     1117           ENDIF !end of ice case 
     1118#endif 
     1119  
     1120        ENDDO !end of loop on the segment 
     1121 
     1122     ENDIF   !end of sec%nb_point =0 case 
     1123     ! 
     1124  END SUBROUTINE transport_h 
     1125  
    8111126  SUBROUTINE dia_dct_sum(sec,jsec)  
    8121127     !!-------------------------------------------------------------  
     
    8901205              SELECT CASE( sec%direction(jseg) )  
    8911206              CASE(0,1)  
    892                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
    893                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    894                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
    895                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
     1207                 ztn   = interp(k%I,k%J,jk,'V',0 )  
     1208                 zsn   = interp(k%I,k%J,jk,'V',1 )  
     1209                 zrhop = interp(k%I,k%J,jk,'V',2)  
     1210                 zrhoi = interp(k%I,k%J,jk,'V',3)  
    8961211 
    8971212              CASE(2,3)  
    898                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
    899                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    900                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
    901                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
     1213                 ztn   = interp(k%I,k%J,jk,'U',0 )  
     1214                 zsn   = interp(k%I,k%J,jk,'U',1 )  
     1215                 zrhop = interp(k%I,k%J,jk,'U',2 )  
     1216                 zrhoi = interp(k%I,k%J,jk,'U',3 )  
    9021217                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    9031218              END SELECT  
     
    9571272                          sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)  
    9581273                       ENDIF  
     1274 
     1275                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1276                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk) 
     1277                       ELSE 
     1278                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk) 
     1279                       ENDIF 
     1280 
     1281                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1282                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk) 
     1283                       ELSE 
     1284                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk) 
     1285                       ENDIF 
    9591286  
    9601287                    ELSE  
     
    9631290                       sec%transport( 5,jclass) = 0._wp  
    9641291                       sec%transport( 6,jclass) = 0._wp  
     1292                       sec%transport( 7,jclass) = 0._wp 
     1293                       sec%transport( 8,jclass) = 0._wp 
     1294                       sec%transport( 9,jclass) = 0._wp 
     1295                       sec%transport(10,jclass) = 0._wp 
    9651296                    ENDIF  
    9661297  
     
    9771308  
    9781309              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN  
    979                  sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6  
     1310                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)*1.E-6  
    9801311              ELSE  
    981                  sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6  
     1312                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)*1.E-6  
    9821313              ENDIF  
    9831314  
    9841315              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN  
    985                  sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6  
     1316                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)*1.E-6  
    9861317              ELSE  
    987                  sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6  
     1318                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)*1.E-6  
    9881319              ENDIF  
    9891320  
     
    9951326     ELSE  !if sec%nb_point =0  
    9961327        sec%transport(1:2,:)=0.  
    997         IF (sec%llstrpond) sec%transport(3:6,:)=0.  
    998         IF (sec%ll_ice_section) sec%transport(7:10,:)=0.  
     1328        IF (sec%llstrpond) sec%transport(3:10,:)=0.  
     1329        IF (sec%ll_ice_section) sec%transport(11:14,:)=0.  
    9991330     ENDIF !end of sec%nb_point =0 case  
    10001331  
    10011332  END SUBROUTINE dia_dct_sum  
     1333 
     1334  SUBROUTINE dia_dct_sum_h(sec,jsec) 
     1335     !!------------------------------------------------------------- 
     1336     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step 
     1337     !! 
     1338     !! Purpose: Average the transport over nn_dctwri time steps  
     1339     !! and sum over the density/salinity/temperature/depth classes 
     1340     !! 
     1341     !! Method:  
     1342     !!           Sum over relevant grid cells to obtain values 
     1343     !!           for each 
     1344     !!              There are several loops:                  
     1345     !!              loop on the segment between 2 nodes 
     1346     !!              loop on the level jk 
     1347     !!              loop on the density/temperature/salinity/level classes 
     1348     !!              test on the density/temperature/salinity/level 
     1349     !! 
     1350     !!  ** Method  :Transport through a given section is equal to the sum of transports 
     1351     !!              computed on each proc. 
     1352     !!              On each proc,transport is equal to the sum of transport computed through 
     1353     !!              segments linking each point of sec%listPoint  with the next one.    
     1354     !! 
     1355     !!------------------------------------------------------------- 
     1356     !! * arguments 
     1357     TYPE(SECTION),INTENT(INOUT) :: sec 
     1358     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     1359 
     1360     TYPE(POINT_SECTION) :: k 
     1361     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes  
     1362     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     1363     !!------------------------------------------------------------- 
     1364 
     1365     !! Sum the relevant segments to obtain values for each class 
     1366     IF(sec%nb_point .NE. 0)THEN    
     1367 
     1368        !--------------------------------------! 
     1369        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     1370        !--------------------------------------! 
     1371        DO jseg=1,MAX(sec%nb_point-1,0) 
     1372            
     1373           !------------------------------------------------------------------------------------------- 
     1374           ! Select the appropriate coordinate for computing the velocity of the segment 
     1375           ! 
     1376           !                      CASE(0)                                    Case (2) 
     1377           !                      -------                                    -------- 
     1378           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     1379           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     1380           !                                                                            | 
     1381           !                                                                            | 
     1382           !                                                                            | 
     1383           !                      Case (3)                                            U(i,j) 
     1384           !                      --------                                              | 
     1385           !                                                                            | 
     1386           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1387           !                        |                                                   | 
     1388           !                        |                                                   | 
     1389           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1390           !                        |                                             
     1391           !                        |                                             
     1392           !                     U(i,j+1)                                             
     1393           !                        |                                       Case(1)      
     1394           !                        |                                       ------       
     1395           !                        |                                             
     1396           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1397           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1398           ! listPoint(jseg)     F(i,j) 
     1399           !  
     1400           !------------------------------------------------------------------------------------------- 
     1401 
     1402           SELECT CASE( sec%direction(jseg) ) 
     1403           CASE(0)  ;   k = sec%listPoint(jseg) 
     1404           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1405           CASE(2)  ;   k = sec%listPoint(jseg) 
     1406           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1407           END SELECT 
     1408 
     1409           !---------------------------| 
     1410           !     LOOP ON THE LEVEL     | 
     1411           !---------------------------| 
     1412           !Sum of the transport on the vertical  
     1413           DO jk=1,mbathy(k%I,k%J) 
     1414 
     1415              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 
     1416              SELECT CASE( sec%direction(jseg) ) 
     1417              CASE(0,1) 
     1418                 ztn   = interp(k%I,k%J,jk,'V',0 ) 
     1419                 zsn   = interp(k%I,k%J,jk,'V',1 ) 
     1420                 zrhop = interp(k%I,k%J,jk,'V',2 ) 
     1421                 zrhoi = interp(k%I,k%J,jk,'V',3 ) 
     1422                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1423              CASE(2,3) 
     1424                 ztn   = interp(k%I,k%J,jk,'U',0 ) 
     1425                 zsn   = interp(k%I,k%J,jk,'U',1 ) 
     1426                 zrhop = interp(k%I,k%J,jk,'U',2 ) 
     1427                 zrhoi = interp(k%I,k%J,jk,'U',3 ) 
     1428                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1429              END SELECT 
     1430 
     1431              zfsdep= fsdept(k%I,k%J,jk) 
     1432  
     1433              !------------------------------- 
     1434              !  LOOP ON THE DENSITY CLASSES | 
     1435              !------------------------------- 
     1436              !The computation is made for each density/heat/salt/... class 
     1437              DO jclass=1,MAX(1,sec%nb_class-1) 
     1438 
     1439                 !----------------------------------------------! 
     1440                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
     1441                 !----------------------------------------------! 
     1442  
     1443                 IF ( (                                                    & 
     1444                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
     1445                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
     1446                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
     1447 
     1448                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
     1449                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
     1450                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
     1451 
     1452                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
     1453                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
     1454                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
     1455 
     1456                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
     1457                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
     1458                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
     1459 
     1460                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
     1461                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
     1462                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
     1463                                                                   ))   THEN 
     1464 
     1465                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 
     1466                    !---------------------------------------------------------------------------- 
     1467                    IF (transports_3d_h(1,jsec,jseg,jk) .GE. 0.0) THEN  
     1468                       sec%transport_h(1,jclass) = sec%transport_h(1,jclass)+transports_3d_h(1,jsec,jseg,jk) 
     1469                    ELSE 
     1470                       sec%transport_h(2,jclass) = sec%transport_h(2,jclass)+transports_3d_h(1,jsec,jseg,jk) 
     1471                    ENDIF 
     1472                    IF( sec%llstrpond )THEN 
     1473 
     1474                       IF ( transports_3d_h(2,jsec,jseg,jk) .GE. 0.0 ) THEN  
     1475                          sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)  
     1476                       ELSE  
     1477                          sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)  
     1478                       ENDIF  
     1479  
     1480                       IF ( transports_3d_h(3,jsec,jseg,jk) .GE. 0.0 ) THEN  
     1481                          sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)  
     1482                       ELSE  
     1483                          sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)  
     1484                       ENDIF  
     1485 
     1486                       IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1487                          sec%transport_h(7,jclass) = sec%transport_h(7,jclass)+transports_3d_h(4,jsec,jseg,jk) 
     1488                       ELSE 
     1489                          sec%transport_h(8,jclass) = sec%transport_h(8,jclass)+transports_3d_h(4,jsec,jseg,jk) 
     1490                       ENDIF 
     1491 
     1492                       IF ( transports_3d_h(5,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1493                          sec%transport_h( 9,jclass) = sec%transport_h( 9,jclass)+transports_3d_h(5,jsec,jseg,jk) 
     1494                       ELSE 
     1495                          sec%transport_h(10,jclass) = sec%transport_h(10,jclass)+transports_3d_h(5,jsec,jseg,jk) 
     1496                       ENDIF 
     1497 
     1498                    ELSE 
     1499                       sec%transport_h( 3,jclass) = 0._wp 
     1500                       sec%transport_h( 4,jclass) = 0._wp 
     1501                       sec%transport_h( 5,jclass) = 0._wp 
     1502                       sec%transport_h( 6,jclass) = 0._wp 
     1503                       sec%transport_h( 7,jclass) = 0._wp 
     1504                       sec%transport_h( 8,jclass) = 0._wp 
     1505                       sec%transport_h( 9,jclass) = 0._wp 
     1506                       sec%transport_h(10,jclass) = 0._wp 
     1507                    ENDIF 
     1508 
     1509                 ENDIF ! end of test if point is in class 
     1510    
     1511              ENDDO ! end of loop on the classes 
     1512 
     1513           ENDDO ! loop over jk 
     1514 
     1515#if defined key_lim2 || defined key_lim3 
     1516 
     1517           !ICE CASE     
     1518           IF( sec%ll_ice_section )THEN 
     1519 
     1520              IF ( transports_2d_h(1,jsec,jseg) .GE. 0.0 ) THEN 
     1521                 sec%transport_h(11,1) = sec%transport_h(11,1)+transports_2d_h(1,jsec,jseg) 
     1522              ELSE 
     1523                 sec%transport_h(12,1) = sec%transport_h(12,1)+transports_2d_h(1,jsec,jseg) 
     1524              ENDIF 
     1525 
     1526              IF ( transports_2d_h(3,jsec,jseg) .GE. 0.0 ) THEN 
     1527                 sec%transport_h(13,1) = sec%transport_h(13,1)+transports_2d_h(2,jsec,jseg) 
     1528              ELSE 
     1529                 sec%transport_h(14,1) = sec%transport_h(14,1)+transports_2d_h(2,jsec,jseg) 
     1530              ENDIF 
     1531 
     1532           ENDIF !end of ice case 
     1533#endif 
     1534  
     1535        ENDDO !end of loop on the segment 
     1536 
     1537     ELSE  !if sec%nb_point =0 
     1538        sec%transport_h(1:2,:)=0. 
     1539        IF (sec%llstrpond) sec%transport_h(3:10,:)=0. 
     1540        IF (sec%ll_ice_section) sec%transport_h( 11:14,:)=0. 
     1541     ENDIF !end of sec%nb_point =0 case 
     1542 
     1543  END SUBROUTINE dia_dct_sum_h 
    10021544   
     1545  SUBROUTINE dia_dct_wri_NOOS(kt,ksec,sec) 
     1546     !!------------------------------------------------------------- 
     1547     !! Write transport output in numdct using NOOS formatting  
     1548     !!  
     1549     !! Purpose: Write  transports in ascii files 
     1550     !!  
     1551     !! Method: 
     1552     !!        1. Write volume transports in "volume_transport" 
     1553     !!           Unit: Sv : area * Velocity / 1.e6  
     1554     !!  
     1555     !!        2. Write heat transports in "heat_transport" 
     1556     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15 
     1557     !!  
     1558     !!        3. Write salt transports in "salt_transport" 
     1559     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6 
     1560     !! 
     1561     !!-------------------------------------------------------------  
     1562     !!arguments 
     1563     INTEGER, INTENT(IN)          :: kt          ! time-step 
     1564     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write    
     1565     INTEGER ,INTENT(IN)          :: ksec        ! section number 
     1566 
     1567     !!local declarations 
     1568     INTEGER               :: jclass,ji             ! Dummy loop 
     1569     CHARACTER(len=2)      :: classe             ! Classname  
     1570     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
     1571     REAL(wp)              :: zslope             ! section's slope coeff 
     1572     ! 
     1573     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
     1574     !!-------------------------------------------------------------  
     1575     CALL wrk_alloc(nb_type_class , zsumclasses )   
     1576 
     1577     zsumclasses(:)=0._wp 
     1578     zslope = sec%slopeSection        
     1579 
     1580     WRITE(numdct_NOOS,'(I4,a1,I2,a1,I2,a12,i3,a17,i3,a10,a25)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1,'   Name: ',sec%name 
     1581 
     1582     DO jclass=1,MAX(1,sec%nb_class-1) 
     1583        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass) 
     1584     ENDDO 
     1585  
     1586     classe   = 'total   ' 
     1587     zbnd1   = 0._wp 
     1588     zbnd2   = 0._wp 
     1589 
     1590     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1591        WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   & 
     1592                                        -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   & 
     1593                                        -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9) 
     1594     ELSE 
     1595        WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   & 
     1596                                          zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   & 
     1597                                          zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10) 
     1598     ENDIF  
     1599 
     1600     DO jclass=1,MAX(1,sec%nb_class-1) 
     1601    
     1602        classe   = 'N       ' 
     1603        zbnd1   = 0._wp 
     1604        zbnd2   = 0._wp 
     1605 
     1606        !insitu density classes transports 
     1607        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. & 
     1608            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN 
     1609           classe = 'DI       ' 
     1610           zbnd1 = sec%zsigi(jclass) 
     1611           zbnd2 = sec%zsigi(jclass+1) 
     1612        ENDIF 
     1613        !potential density classes transports 
     1614        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. & 
     1615            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN 
     1616           classe = 'DP      ' 
     1617           zbnd1 = sec%zsigp(jclass) 
     1618           zbnd2 = sec%zsigp(jclass+1) 
     1619        ENDIF 
     1620        !depth classes transports 
     1621        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. & 
     1622            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN  
     1623           classe = 'Z       ' 
     1624           zbnd1 = sec%zlay(jclass) 
     1625           zbnd2 = sec%zlay(jclass+1) 
     1626        ENDIF 
     1627        !salinity classes transports 
     1628        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. & 
     1629            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN 
     1630           classe = 'S       ' 
     1631           zbnd1 = sec%zsal(jclass) 
     1632           zbnd2 = sec%zsal(jclass+1)    
     1633        ENDIF 
     1634        !temperature classes transports 
     1635        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. & 
     1636            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN 
     1637           classe = 'T       ' 
     1638           zbnd1 = sec%ztem(jclass) 
     1639           zbnd2 = sec%ztem(jclass+1) 
     1640        ENDIF 
     1641                   
     1642        !write volume transport per class 
     1643        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1644           WRITE(numdct_NOOS,'(9e12.4E2)') -(sec%transport( 1,jclass)+sec%transport( 2,jclass)),-sec%transport( 2,jclass),-sec%transport( 1,jclass), & 
     1645                                           -(sec%transport( 7,jclass)+sec%transport( 8,jclass)),-sec%transport( 8,jclass),-sec%transport( 7,jclass), & 
     1646                                           -(sec%transport( 9,jclass)+sec%transport(10,jclass)),-sec%transport(10,jclass),-sec%transport( 9,jclass) 
     1647        ELSE 
     1648           WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), & 
     1649                                             sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), & 
     1650                                             sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass) 
     1651        ENDIF 
     1652 
     1653     ENDDO 
     1654 
     1655     CALL wrk_dealloc(nb_type_class , zsumclasses )   
     1656 
     1657  END SUBROUTINE dia_dct_wri_NOOS 
     1658 
     1659  SUBROUTINE dia_dct_wri_NOOS_h(hr,ksec,sec) 
     1660     !!------------------------------------------------------------- 
     1661     !! As routine dia_dct_wri_NOOS but for hourly output files 
     1662     !! 
     1663     !! Write transport output in numdct using NOOS formatting  
     1664     !!  
     1665     !! Purpose: Write  transports in ascii files 
     1666     !!  
     1667     !! Method: 
     1668     !!        1. Write volume transports in "volume_transport" 
     1669     !!           Unit: Sv : area * Velocity / 1.e6  
     1670     !! 
     1671     !!-------------------------------------------------------------  
     1672     !!arguments 
     1673     INTEGER, INTENT(IN)          :: hr          ! hour 
     1674     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write    
     1675     INTEGER ,INTENT(IN)          :: ksec        ! section number 
     1676 
     1677     !!local declarations 
     1678     INTEGER               :: jclass,jhr            ! Dummy loop 
     1679     CHARACTER(len=2)      :: classe             ! Classname  
     1680     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
     1681     REAL(wp)              :: zslope             ! section's slope coeff 
     1682     ! 
     1683     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
     1684     !!-------------------------------------------------------------  
     1685 
     1686     CALL wrk_alloc(nb_type_class , zsumclasses )  
     1687 
     1688     zsumclasses(:)=0._wp 
     1689     zslope = sec%slopeSection        
     1690 
     1691     DO jclass=1,MAX(1,sec%nb_class-1) 
     1692        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport_h(1:nb_type_class,jclass) 
     1693     ENDDO 
     1694  
     1695     !write volume transport per class 
     1696     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1697        z_hr_output(ksec,hr,1)=-(zsumclasses(1)+zsumclasses(2)) 
     1698     ELSE 
     1699        z_hr_output(ksec,hr,1)= (zsumclasses(1)+zsumclasses(2)) 
     1700     ENDIF 
     1701 
     1702     DO jclass=1,MAX(1,sec%nb_class-1) 
     1703        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1704           z_hr_output(ksec,hr,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 
     1705        ELSE 
     1706           z_hr_output(ksec,hr,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 
     1707        ENDIF 
     1708     ENDDO 
     1709 
     1710     IF ( hr .eq. 48._wp ) THEN 
     1711        WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1 
     1712        DO jhr=25,48 
     1713           WRITE(numdct_NOOS_h,'(11F12.1)')  z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) ) 
     1714        ENDDO 
     1715     ENDIF 
     1716 
     1717     CALL wrk_dealloc(nb_type_class , zsumclasses ) 
     1718 
     1719  END SUBROUTINE dia_dct_wri_NOOS_h 
     1720 
    10031721  SUBROUTINE dia_dct_wri(kt,ksec,sec) 
    10041722     !!------------------------------------------------------------- 
     
    10921810           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  & 
    10931811                              jclass,classe,zbnd1,zbnd2,& 
    1094                               sec%transport(3,jclass)*1.e-15,sec%transport(4,jclass)*1.e-15, & 
    1095                               ( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15 
     1812                              sec%transport(7,jclass)*1.e-15,sec%transport(8,jclass)*1.e-15, & 
     1813                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1.e-15 
    10961814           !write salt transport per class 
    10971815           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  & 
    10981816                              jclass,classe,zbnd1,zbnd2,& 
    1099                               sec%transport(5,jclass)*1.e-9,sec%transport(6,jclass)*1.e-9,& 
    1100                               (sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9 
     1817                              sec%transport(9,jclass)*1.e-9,sec%transport(10,jclass)*1.e-9,& 
     1818                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1.e-9 
    11011819        ENDIF 
    11021820 
     
    11171835        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, & 
    11181836                           jclass,"total",zbnd1,zbnd2,& 
    1119                            zsumclasses(3)*1.e-15,zsumclasses(4)*1.e-15,& 
    1120                            (zsumclasses(3)+zsumclasses(4) )*1.e-15 
     1837                           zsumclasses(7)*1.e-15,zsumclasses(8)*1.e-15,& 
     1838                           (zsumclasses(7)+zsumclasses(8) )*1.e-15 
    11211839        !write total salt transport 
    11221840        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, & 
    11231841                           jclass,"total",zbnd1,zbnd2,& 
    1124                            zsumclasses(5)*1.e-9,zsumclasses(6)*1.e-9,& 
    1125                            (zsumclasses(5)+zsumclasses(6))*1.e-9 
     1842                           zsumclasses(9)*1.e-9,zsumclasses(10)*1.e-9,& 
     1843                           (zsumclasses(9)+zsumclasses(10))*1.e-9 
    11261844     ENDIF 
    11271845 
     
    11311849        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 
    11321850                              jclass,"ice_vol",zbnd1,zbnd2,& 
    1133                               sec%transport(7,1),sec%transport(8,1),& 
    1134                               sec%transport(7,1)+sec%transport(8,1) 
     1851                              sec%transport(11,1),sec%transport(12,1),& 
     1852                              sec%transport(11,1)+sec%transport(12,1) 
    11351853        !write total ice surface transport 
    11361854        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 
    11371855                              jclass,"ice_surf",zbnd1,zbnd2,& 
    1138                               sec%transport(9,1),sec%transport(10,1), & 
    1139                               sec%transport(9,1)+sec%transport(10,1)  
     1856                              sec%transport(13,1),sec%transport(14,1), & 
     1857                              sec%transport(13,1)+sec%transport(14,1)  
    11401858     ENDIF 
    11411859                                               
     
    11461864  END SUBROUTINE dia_dct_wri 
    11471865 
    1148   FUNCTION interp(ki, kj, kk, cd_point, ptab) 
     1866   PURE  FUNCTION interp(ki, kj, kk, cd_point,var)  
    11491867  !!---------------------------------------------------------------------- 
    11501868  !! 
     
    12081926  !*arguments 
    12091927  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point 
     1928  INTEGER, INTENT(IN)                          :: var   !  which variable 
    12101929  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V) 
    1211   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk ) 
    12121930  REAL(wp)                                     :: interp       ! interpolated variable  
    12131931 
     
    12201938  !!---------------------------------------------------------------------- 
    12211939 
     1940  
     1941 
    12221942  IF( cd_point=='U' )THEN  
    12231943     ii1 = ki    ; ij1 = kj  
     
    12501970   
    12511971     ! result 
    1252      interp = zmsk * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )    
    1253  
     1972           SELECT CASE( var ) 
     1973              CASE(0)  ;    interp = zmsk * ( zwgt2 *  tsn(ii1,ij1,kk,jp_tem) + zwgt1 *  tsn(ii1,ij1,kk,jp_tem) ) / ( zwgt2 + zwgt1 ) 
     1974              CASE(1)  ;    interp = zmsk * ( zwgt2 *  tsn(ii1,ij1,kk,jp_sal) + zwgt1 *  tsn(ii1,ij1,kk,jp_sal) ) / ( zwgt2 + zwgt1 ) 
     1975              CASE(2)  ;    interp = zmsk * ( zwgt2 *  rhop(ii1,ij1,kk) + zwgt1 *  rhop(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 ) 
     1976              CASE(3)  ;    interp = zmsk * ( zwgt2 *  (rhd(ii1,ij1,kk)*rau0+rau0) + zwgt1 *  (rhd(ii1,ij1,kk)*rau0+rau0) ) / ( zwgt2 + zwgt1 ) 
     1977           END SELECT 
    12541978 
    12551979  ELSE       ! full step or partial step case  
     
    12731997        IF( ze3t >= 0. )THEN  
    12741998           ! zbis 
    1275            zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )  
     1999           SELECT CASE( var ) 
     2000           CASE(0)   
     2001                     zbis   =  tsn(ii2,ij2,kk,jp_tem) + zwgt1 *  (tsn(ii2,ij2,kk-1,jp_tem)-tsn(ii2,ij2,kk,jp_tem)   ) 
     2002                     interp =  zmsk * ( zet2 * tsn(ii1,ij1,kk,jp_tem) + zet1 * zbis )/( zet1 + zet2 ) 
     2003           CASE(1)   
     2004                     zbis   =  tsn(ii2,ij2,kk,jp_sal) + zwgt1 *  (tsn(ii2,ij2,kk-1,jp_sal)-tsn(ii2,ij2,kk,jp_sal)   ) 
     2005                     interp =  zmsk * ( zet2 * tsn(ii1,ij1,kk,jp_sal) + zet1 * zbis )/( zet1 + zet2 ) 
     2006           CASE(2)   
     2007                     zbis   =  rhop(ii2,ij2,kk) + zwgt1 *  (rhop(ii2,ij2,kk-1)-rhop(ii2,ij2,kk)   ) 
     2008                     interp =  zmsk * ( zet2 * rhop(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 
     2009           CASE(3)   
     2010                     zbis   =  (rhd(ii2,ij2,kk)*rau0+rau0) + zwgt1 *  ( (rhd(ii2,ij2,kk-1)*rau0+rau0)-(rhd(ii2,ij2,kk)*rau0+rau0)   ) 
     2011                     interp =  zmsk * ( zet2 * (rhd(ii1,ij1,kk)*rau0+rau0) + zet1 * zbis )/( zet1 + zet2 ) 
     2012           END SELECT 
    12762013           ! result 
    1277             interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 
    12782014        ELSE 
    12792015           ! zbis 
    1280            zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) ) 
    1281            ! result 
    1282            interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2016           SELECT CASE( var ) 
     2017           CASE(0)   
     2018                 zbis   = tsn(ii1,ij1,kk,jp_tem) + zwgt2 * ( tsn(ii1,ij1,kk-1,jp_tem) - tsn(ii1,ij2,kk,jp_tem) ) 
     2019                 interp = zmsk * ( zet2 * zbis + zet1 * tsn(ii2,ij2,kk,jp_tem) )/( zet1 + zet2 ) 
     2020           CASE(1)   
     2021                 zbis   = tsn(ii1,ij1,kk,jp_sal) + zwgt2 * ( tsn(ii1,ij1,kk-1,jp_sal) - tsn(ii1,ij2,kk,jp_sal) ) 
     2022                 interp = zmsk * ( zet2 * zbis + zet1 * tsn(ii2,ij2,kk,jp_sal) )/( zet1 + zet2 ) 
     2023           CASE(2)   
     2024                 zbis   = rhop(ii1,ij1,kk) + zwgt2 * ( rhop(ii1,ij1,kk-1) - rhop(ii1,ij2,kk) ) 
     2025                 interp = zmsk * ( zet2 * zbis + zet1 * rhop(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2026           CASE(3)   
     2027                 zbis   = (rhd(ii1,ij1,kk)*rau0+rau0) + zwgt2 * ( (rhd(ii1,ij1,kk-1)*rau0+rau0) - (rhd(ii1,ij2,kk)*rau0+rau0) ) 
     2028                 interp = zmsk * ( zet2 * zbis + zet1 * (rhd(ii2,ij2,kk)*rau0+rau0) )/( zet1 + zet2 ) 
     2029           END SELECT 
    12832030        ENDIF     
    12842031 
    12852032     ELSE 
    1286         interp = zmsk * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2033        SELECT CASE( var ) 
     2034        CASE(0)   
     2035           interp = zmsk * (  zet2 * tsn(ii1,ij1,kk,jp_tem) + zet1 * tsn(ii2,ij2,kk,jp_tem) )/( zet1 + zet2 ) 
     2036        CASE(1)   
     2037           interp = zmsk * (  zet2 * tsn(ii1,ij1,kk,jp_sal) + zet1 * tsn(ii2,ij2,kk,jp_sal) )/( zet1 + zet2 ) 
     2038        CASE(2)   
     2039           interp = zmsk * (  zet2 * rhop(ii1,ij1,kk) + zet1 * rhop(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2040        CASE(3)   
     2041           interp = zmsk * (  zet2 * (rhd(ii1,ij1,kk)*rau0+rau0) + zet1 * (rhd(ii2,ij2,kk)*rau0+rau0) )/( zet1 + zet2 ) 
     2042        END SELECT 
    12872043     ENDIF 
    12882044 
    12892045  ENDIF 
    1290  
    12912046 
    12922047  END FUNCTION interp 
Note: See TracChangeset for help on using the changeset viewer.