source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 10237

Last change on this file since 10237 was 10237, checked in by kingr, 2 years ago

Merged in NOOS transect code changes from AMM15_v3_6_STABLE_package_withNOOS.

  • Property svn:executable set to *
File size: 97.8 KB
Line 
1MODULE 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  !!----------------------------------------------------------------------
18#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
40#if defined key_lim2
41  USE ice_2
42#endif
43#if defined key_lim3
44  USE ice
45#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  PRIVATE  dia_dct_wri_NOOS
62
63#include "domzgr_substitute.h90"
64
65  !! * Shared module variables
66  LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
67
68  !! * Module variables
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
74   
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
81  INTEGER            :: nb_sec 
82
83  TYPE POINT_SECTION
84     INTEGER :: I,J
85  END TYPE POINT_SECTION
86
87  TYPE COORD_SECTION
88     REAL(wp) :: lon,lat
89  END TYPE COORD_SECTION
90
91  TYPE SECTION
92     CHARACTER(len=60)                            :: name              ! name of the sec
93     LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
94                                                                       ! heat transports
95     LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
96     LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
97     TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
98     INTEGER                                      :: nb_class          ! number of boundaries for density classes
99     INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
100     CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
101     REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
102                                                     zsigp           ,&! potential density classes    (99 if you don't want)
103                                                     zsal            ,&! salinity classes   (99 if you don't want)
104                                                     ztem            ,&! temperature classes(99 if you don't want)
105                                                     zlay              ! level classes      (99 if you don't want)
106     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
108     REAL(wp)                                         :: slopeSection  ! slope of the section
109     INTEGER                                          :: nb_point      ! number of points in the section
110     TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
111  END TYPE SECTION
112
113  TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
114 
115  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
116  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
120
121   !! $Id$
122CONTAINS
123
124 
125  INTEGER FUNCTION diadct_alloc() 
126     !!----------------------------------------------------------------------
127     !!                   ***  FUNCTION diadct_alloc  ***
128     !!----------------------------------------------------------------------
129     INTEGER :: ierr(5) 
130     !!----------------------------------------------------------------------
131
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
137
138     diadct_alloc = MAXVAL( ierr ) 
139     IF( diadct_alloc /= 0 )   CALL ctl_warn('diadct_alloc: failed to allocate arrays') 
140 
141  END FUNCTION diadct_alloc 
142
143  SUBROUTINE dia_dct_init
144     !!---------------------------------------------------------------------
145     !!               ***  ROUTINE diadct  *** 
146     !!
147     !!  ** Purpose: Read the namelist parameters
148     !!              Open output files
149     !!
150     !!---------------------------------------------------------------------
151     NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
152     INTEGER  ::   ios                 ! Local integer output status for namelist read
153
154     IF( nn_timing == 1 )   CALL timing_start('dia_dct_init')
155
156     REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections
157     READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
158901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp )
159
160     REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
161     READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
162902  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp )
163     IF(lwm) WRITE ( numond, namdct )
164
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
172     IF( lwp ) THEN
173        WRITE(numout,*) " "
174        WRITE(numout,*) "diadct_init: compute transports through sections "
175        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
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
185
186        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
187                                            WRITE(numout,*)"       Debug section number: ", nn_secdebug 
188        ELSE IF ( nn_secdebug ==  0 )THEN ; WRITE(numout,*)"       No section to debug"
189        ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)"       Debug all sections"
190        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
191        ENDIF
192
193        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
194          &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
195
196     ENDIF
197
198     !Read section_ijglobal.diadct
199     CALL readsec
200
201     !open output file
202     IF( lwm ) THEN
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
211     ENDIF
212
213     ! Initialise arrays to zero
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
219
220     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init')
221     !
222  END SUBROUTINE dia_dct_init
223 
224 
225  SUBROUTINE dia_dct(kt)
226     !!---------------------------------------------------------------------
227     !!               ***  ROUTINE diadct  *** 
228     !!
229     !!  Purpose :: Compute section transports and write it in numdct files
230     !!   
231     !!  Method  :: All arrays initialised to zero in dct_init
232     !!             Each nn_dct time step call subroutine 'transports' for
233     !!               each section to sum the transports over each grid cell.
234     !!             Each nn_dctwri time step:
235     !!               Divide the arrays by the number of summations to gain
236     !!               an average value
237     !!               Call dia_dct_sum to sum relevant grid boxes to obtain
238     !!               totals for each class (density, depth, temp or sal)
239     !!               Call dia_dct_wri to write the transports into file
240     !!               Reinitialise all relevant arrays to zero
241     !!---------------------------------------------------------------------
242     !! * Arguments
243     INTEGER,INTENT(IN)        ::kt
244
245     !! * Local variables
246     INTEGER             :: jsec,            &! loop on sections
247                            itotal            ! nb_sec_max*nb_type_class*nb_class_max
248     LOGICAL             :: lldebug =.FALSE.  ! debug a section 
249     
250     INTEGER , DIMENSION(1)             :: ish   ! tmp array for mpp_sum
251     INTEGER , DIMENSION(3)             :: ish2  !   "
252     REAL(wp), POINTER, DIMENSION(:)    :: zwork !   " 
253     REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   "
254
255     !!---------------------------------------------------------------------   
256     IF( nn_timing == 1 )   CALL timing_start('dia_dct')
257
258     IF( lk_mpp )THEN
259        itotal = nb_sec_max*nb_type_class*nb_class_max
260        CALL wrk_alloc( itotal                                , zwork ) 
261        CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum  )
262     ENDIF   
263     lldebug=.TRUE.
264     ! Initialise arrays
265     zwork(:) = 0.0 
266     zsum(:,:,:) = 0.0
267
268     IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
269         WRITE(numout,*) " "
270         WRITE(numout,*) "diadct: compute transport"
271         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
272         WRITE(numout,*) "nb_sec = ",nb_sec
273     ENDIF
274
275 
276     ! Compute transport and write only at nn_dctwri
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
279
280        DO jsec=1,nb_sec
281
282           !debug this section computing ?
283           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. 
284
285           !Compute transport through section 
286           CALL transport(secs(jsec),lldebug,jsec) 
287
288        ENDDO
289             
290        IF( MOD(kt,nn_dctwri)==0 )THEN
291
292           IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: average transports and write at kt = ",kt         
293 
294           !! divide arrays by nn_dctwri/nn_dct to obtain average
295           transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 
296           transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct) 
297 
298           ! Sum over each class
299           DO jsec=1,nb_sec 
300              CALL dia_dct_sum(secs(jsec),jsec) 
301           ENDDO 
302
303           !Sum on all procs
304           IF( lk_mpp )THEN
305              zsum(:,:,:)=0.0_wp
306              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
307              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
308              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
309              zwork(:)= RESHAPE(zsum(:,:,:), ish )
310              CALL mpp_sum(zwork, ish(1))
311              zsum(:,:,:)= RESHAPE(zwork,ish2)
312              DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
313           ENDIF
314
315           !Write the transport
316           DO jsec=1,nb_sec
317
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
320           
321              !nullify transports values after writing
322              transports_3d(:,jsec,:,:)=0.
323              transports_2d(:,jsec,:  )=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)
326
327           ENDDO
328
329        ENDIF
330
331     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   
385
386     IF( lk_mpp )THEN
387        itotal = nb_sec_max*nb_type_class*nb_class_max
388        CALL wrk_dealloc( itotal                                , zwork ) 
389        CALL wrk_dealloc( nb_sec_max,nb_type_class,nb_class_max , zsum  )
390     ENDIF   
391
392     IF( nn_timing == 1 )   CALL timing_stop('dia_dct')
393     !
394  END SUBROUTINE dia_dct
395
396  SUBROUTINE readsec 
397     !!---------------------------------------------------------------------
398     !!               ***  ROUTINE readsec  ***
399     !!
400     !!  ** Purpose:
401     !!            Read a binary file(section_ijglobal.diadct)
402     !!            generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
403     !!
404     !!
405     !!---------------------------------------------------------------------
406     !! * Local variables
407     INTEGER :: iptglo , iptloc                               ! Global and local number of points for a section
408     INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2  ! temporary  integer
409     INTEGER :: jsec, jpt                                     ! dummy loop indices
410
411     INTEGER, DIMENSION(2) :: icoord 
412     CHARACTER(len=160)    :: clname                          !filename
413     CHARACTER(len=200)    :: cltmp
414     CHARACTER(len=200)    :: clformat                        !automatic format
415     TYPE(POINT_SECTION),DIMENSION(nb_point_max)  ::coordtemp !contains listpoints coordinates
416                                                              !read in the file
417     INTEGER, POINTER, DIMENSION(:) :: directemp              !contains listpoints directions
418                                                              !read in the files
419     LOGICAL :: llbon                                       ,&!local logical
420                lldebug                                       !debug the section
421     !!-------------------------------------------------------------------------------------
422     CALL wrk_alloc( nb_point_max, directemp )
423
424     !open input file
425     !---------------
426     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
427 
428     !---------------
429     !Read input file
430     !---------------
431     
432     DO jsec=1,nb_sec_max      !loop on the nb_sec sections
433
434        IF ( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) ) &
435           & WRITE(numout,*)'debuging for section number: ',jsec 
436
437        !initialization
438        !---------------
439        secs(jsec)%name=''
440        secs(jsec)%llstrpond    = .FALSE.  ; secs(jsec)%ll_ice_section = .FALSE.
441        secs(jsec)%ll_date_line = .FALSE.  ; secs(jsec)%nb_class       = 0
442        secs(jsec)%zsigi        = 99._wp   ; secs(jsec)%zsigp          = 99._wp
443        secs(jsec)%zsal         = 99._wp   ; secs(jsec)%ztem           = 99._wp
444        secs(jsec)%zlay         = 99._wp         
445        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0
446        secs(jsec)%transport_h  =  0._wp   ; secs(jsec)%nb_point       = 0
447
448        !read section's number / name / computing choices / classes / slopeSection / points number
449        !-----------------------------------------------------------------------------------------
450        READ(numdct_in,iostat=iost)isec
451        IF (iost .NE. 0 )EXIT !end of file
452        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
453        IF( jsec .NE. isec )CALL ctl_stop( cltmp )
454
455        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec 
456
457        READ(numdct_in)secs(jsec)%name
458        READ(numdct_in)secs(jsec)%llstrpond
459        READ(numdct_in)secs(jsec)%ll_ice_section
460        READ(numdct_in)secs(jsec)%ll_date_line
461        READ(numdct_in)secs(jsec)%coordSec
462        READ(numdct_in)secs(jsec)%nb_class
463        READ(numdct_in)secs(jsec)%zsigi
464        READ(numdct_in)secs(jsec)%zsigp
465        READ(numdct_in)secs(jsec)%zsal
466        READ(numdct_in)secs(jsec)%ztem
467        READ(numdct_in)secs(jsec)%zlay
468        READ(numdct_in)secs(jsec)%slopeSection
469        READ(numdct_in)iptglo
470
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
477        !debug
478        !-----
479
480        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
481         
482            WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))' 
483
484            WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name)
485            WRITE(numout,*)       "      Compute heat and salt transport ? ",secs(jsec)%llstrpond
486            WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section
487            WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line
488            WRITE(numout,*)       "      Slope section :                   ",secs(jsec)%slopeSection
489            WRITE(numout,*)       "      Number of points in the section:  ",iptglo
490            WRITE(numout,*)       "      Number of classes                 ",secs(jsec)%nb_class
491            WRITE(numout,clformat)"      Insitu density classes :          ",secs(jsec)%zsigi
492            WRITE(numout,clformat)"      Potential density classes :       ",secs(jsec)%zsigp
493            WRITE(numout,clformat)"      Salinity classes :                ",secs(jsec)%zsal
494            WRITE(numout,clformat)"      Temperature classes :             ",secs(jsec)%ztem
495            WRITE(numout,clformat)"      Depth classes :                   ",secs(jsec)%zlay
496        ENDIF               
497
498        IF( iptglo .NE. 0 )THEN
499             
500           !read points'coordinates and directions
501           !--------------------------------------
502           IF ( ln_NOOS ) THEN
503              WRITE(numout,*) 'Coords and direction read'
504           ENDIF
505
506           coordtemp(:) = POINT_SECTION(0,0) !list of points read
507           directemp(:) = 0                  !value of directions of each points
508           DO jpt=1,iptglo
509              READ(numdct_in)i1,i2
510              coordtemp(jpt)%I = i1 
511              coordtemp(jpt)%J = i2
512           ENDDO
513           READ(numdct_in)directemp(1:iptglo)
514   
515           !debug
516           !-----
517           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
518              WRITE(numout,*)"      List of points in global domain:"
519              DO jpt=1,iptglo
520                 WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt),directemp(jpt)
521              ENDDO                 
522           ENDIF
523
524           !Now each proc selects only points that are in its domain:
525           !--------------------------------------------------------
526           iptloc = 0                    !initialize number of points selected
527           DO jpt=1,iptglo               !loop on listpoint read in the file
528                   
529              iiglo=coordtemp(jpt)%I          ! global coordinates of the point
530              ijglo=coordtemp(jpt)%J          !  "
531
532              IF( iiglo==jpidta .AND. nimpp==1 ) iiglo = 2
533
534              iiloc=iiglo-jpizoom+1-nimpp+1   ! local coordinates of the point
535              ijloc=ijglo-jpjzoom+1-njmpp+1   !  "
536
537              !verify if the point is on the local domain:(1,nlei)*(1,nlej)
538              IF( iiloc .GE. 1 .AND. iiloc .LE. nlei .AND. &
539                  ijloc .GE. 1 .AND. ijloc .LE. nlej       )THEN
540                 iptloc = iptloc + 1                                                 ! count local points
541                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
542                 secs(jsec)%direction(iptloc) = directemp(jpt)                       ! store local direction
543              ENDIF
544
545           ENDDO
546     
547           secs(jsec)%nb_point=iptloc !store number of section's points
548
549           !debug
550           !-----
551           IF(   lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
552              WRITE(numout,*)"      List of points selected by the proc:"
553              DO jpt = 1,iptloc
554                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
555                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
556                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
557              ENDDO
558           ENDIF
559
560              IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
561              DO jpt = 1,iptloc
562                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
563                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
564              ENDDO
565              ENDIF
566
567           !remove redundant points between processors
568           !------------------------------------------
569           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
570           IF( iptloc .NE. 0 )THEN
571              CALL removepoints(secs(jsec),'I','top_list',lldebug)
572              CALL removepoints(secs(jsec),'I','bot_list',lldebug)
573              CALL removepoints(secs(jsec),'J','top_list',lldebug)
574              CALL removepoints(secs(jsec),'J','bot_list',lldebug)
575           ENDIF
576           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
577              DO jpt = 1,secs(jsec)%nb_point
578                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
579                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
580              ENDDO
581           ENDIF
582
583           !debug
584           !-----
585           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
586              WRITE(numout,*)"      List of points after removepoints:"
587              iptloc = secs(jsec)%nb_point
588              DO jpt = 1,iptloc
589                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
590                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
591                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
592                 CALL FLUSH(numout)
593              ENDDO
594           ENDIF
595
596        ELSE  ! iptglo = 0
597           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )&
598              WRITE(numout,*)'   No points for this section.'
599        ENDIF
600
601     ENDDO !end of the loop on jsec
602 
603     nb_sec = jsec-1   !number of section read in the file
604
605     CALL wrk_dealloc( nb_point_max, directemp )
606     !
607  END SUBROUTINE readsec
608
609  SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
610     !!---------------------------------------------------------------------------
611     !!             *** function removepoints
612     !!
613     !!   ** Purpose :: Remove points which are common to 2 procs
614     !!
615     !----------------------------------------------------------------------------
616     !! * arguments
617     TYPE(SECTION),INTENT(INOUT) :: sec
618     CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
619     CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
620     LOGICAL,INTENT(IN)          :: ld_debug                     
621
622     !! * Local variables
623     INTEGER :: iextr         ,& !extremity of listpoint that we verify
624                iind          ,& !coord     of listpoint that we verify
625                itest         ,& !indice value of the side of the domain
626                                 !where points could be redundant
627                isgn          ,& ! isgn= 1 : scan listpoint from start to end
628                                 ! isgn=-1 : scan listpoint from end to start
629                istart,iend      !first and last points selected in listpoint
630     INTEGER :: jpoint           !loop on list points
631     INTEGER, POINTER, DIMENSION(:)   :: idirec !contains temporary sec%direction
632     INTEGER, POINTER, DIMENSION(:,:) :: icoord !contains temporary sec%listpoint
633     !----------------------------------------------------------------------------
634     CALL wrk_alloc(    nb_point_max, idirec )
635     CALL wrk_alloc( 2, nb_point_max, icoord )
636
637     IF( ld_debug )WRITE(numout,*)'      -------------------------'
638     IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
639
640     !iextr=extremity of list_point that we verify
641     IF      ( cdextr=='bot_list' )THEN ; iextr=1            ; isgn=1
642     ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
643     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdextr")
644     ENDIF
645 
646     !which coordinate shall we verify ?
647     IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1
648     ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2
649     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind") 
650     ENDIF
651
652     IF( ld_debug )THEN
653        WRITE(numout,*)'      case: coord/list extr/domain side'
654        WRITE(numout,*)'      ', cdind,' ',cdextr,' ',itest
655        WRITE(numout,*)'      Actual number of points: ',sec%nb_point
656     ENDIF
657
658     icoord(1,1:nb_point_max) = sec%listPoint%I
659     icoord(2,1:nb_point_max) = sec%listPoint%J
660     idirec                   = sec%direction
661     sec%listPoint            = POINT_SECTION(0,0)
662     sec%direction            = 0
663
664     jpoint=iextr+isgn
665     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
666         IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
667         ELSE                                                                               ; EXIT
668         ENDIF
669     ENDDO 
670
671     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
672     ELSE                        ; istart=1        ; iend=jpoint+1
673     ENDIF
674
675     sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
676     sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
677     sec%direction(1:1+iend-istart)   = idirec(istart:iend)
678     sec%nb_point                     = iend-istart+1
679     
680     IF( ld_debug )THEN
681        WRITE(numout,*)'      Number of points after removepoints :',sec%nb_point
682        WRITE(numout,*)'      sec%direction after removepoints :',sec%direction(1:sec%nb_point)
683     ENDIF
684
685     CALL wrk_dealloc(    nb_point_max, idirec )
686     CALL wrk_dealloc( 2, nb_point_max, icoord )
687  END SUBROUTINE removepoints
688
689  SUBROUTINE transport(sec,ld_debug,jsec)
690     !!-------------------------------------------------------------------------------------------
691     !!                     ***  ROUTINE transport  ***
692     !!
693     !!  Purpose ::  Compute the transport for each point in a section
694     !!
695     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
696     !!              Be aware :           
697     !!              One section is a sum of segments
698     !!              One segment is defined by 2 consecutive points in sec%listPoint
699     !!              All points of sec%listPoint are positioned on the F-point of the cell
700     !!
701     !!              There are two loops:                 
702     !!              loop on the segment between 2 nodes
703     !!              loop on the level jk !!
704     !!
705     !!  Output  ::  Arrays containing the volume,density,heat,salt transports for each i
706     !!              point in a section, summed over each nn_dct.
707     !!
708     !!-------------------------------------------------------------------------------------------
709     !! * Arguments
710     TYPE(SECTION),INTENT(INOUT) :: sec
711     LOGICAL      ,INTENT(IN)    :: ld_debug
712     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
713   
714     !! * Local variables
715     INTEGER             :: jk, jseg, jclass,jl,                 &!loop on level/segment/classes/ice categories
716                            isgnu, isgnv                          !
717     REAL(wp)            :: zumid, zvmid,                        &!U/V velocity on a cell segment
718                            zumid_ice, zvmid_ice,                &!U/V ice velocity
719                            zTnorm                                !transport of velocity through one cell's sides
720     REAL(wp)            :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point
721
722     TYPE(POINT_SECTION) :: k
723     !!--------------------------------------------------------
724
725     !!NIALL IF( ld_debug )WRITE(numout,*)'      Compute transport'
726
727     !---------------------------!
728     !  COMPUTE TRANSPORT        !
729     !---------------------------!
730     IF(sec%nb_point .NE. 0)THEN   
731
732        !----------------------------------------------------------------------------------------------------
733        !Compute sign for velocities:
734        !
735        !convention:
736        !   non horizontal section: direction + is toward left hand of section
737        !       horizontal section: direction + is toward north of section
738        !
739        !
740        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
741        !       ----------------      -----------------     ---------------             --------------
742        !
743        !   isgnv=1         direction +     
744        !  ______         _____             ______                                                   
745        !        |           //|            |                  |                         direction +   
746        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
747        !        |_______  //         ______|    \\            | ---\                        |
748        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
749        !               |             |          __\\|         |                   
750        !               |             |     direction +        |                      isgnv=1                                 
751        !                                                     
752        !----------------------------------------------------------------------------------------------------
753        isgnu = 1
754        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
755        ELSE                                ; isgnv =  1
756        ENDIF
757        IF( sec%slopeSection .GE. 9999. )     isgnv =  1
758
759        IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
760
761        !--------------------------------------!
762        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
763        !--------------------------------------!
764        DO jseg=1,MAX(sec%nb_point-1,0)
765             
766           !-------------------------------------------------------------------------------------------
767           ! Select the appropriate coordinate for computing the velocity of the segment
768           !
769           !                      CASE(0)                                    Case (2)
770           !                      -------                                    --------
771           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
772           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
773           !                                                                            |
774           !                                                                            |
775           !                                                                            |
776           !                      Case (3)                                            U(i,j)
777           !                      --------                                              |
778           !                                                                            |
779           !  listPoint(jseg+1) F(i,j+1)                                                |
780           !                        |                                                   |
781           !                        |                                                   |
782           !                        |                                 listPoint(jseg+1) F(i,j-1)
783           !                        |                                           
784           !                        |                                           
785           !                     U(i,j+1)                                           
786           !                        |                                       Case(1)     
787           !                        |                                       ------     
788           !                        |                                           
789           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
790           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
791           ! listPoint(jseg)     F(i,j)
792           !
793           !-------------------------------------------------------------------------------------------
794
795           SELECT CASE( sec%direction(jseg) )
796           CASE(0)  ;   k = sec%listPoint(jseg)
797           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
798           CASE(2)  ;   k = sec%listPoint(jseg)
799           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
800           END SELECT
801
802           !---------------------------|
803           !     LOOP ON THE LEVEL     |
804           !---------------------------|
805           !Sum of the transport on the vertical 
806           DO jk=1,mbathy(k%I,k%J) 
807 
808              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
809              SELECT CASE( sec%direction(jseg) ) 
810              CASE(0,1) 
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 ) 
815                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
816              CASE(2,3) 
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 ) 
821                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
822              END SELECT
823 
824              zfsdep= fsdept(k%I,k%J,jk) 
825 
826              !compute velocity with the correct direction
827              SELECT CASE( sec%direction(jseg) ) 
828              CASE(0,1)   
829                 zumid=0. 
830                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
831              CASE(2,3) 
832                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
833                 zvmid=0. 
834              END SELECT 
835 
836              !zTnorm=transport through one cell;
837              !velocity* cell's length * cell's thickness
838              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
839                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
840
841#if ! defined key_vvl
842              !add transport due to free surface
843              IF( jk==1 )THEN
844                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
845                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
846              ENDIF 
847#endif
848              !COMPUTE TRANSPORT 
849 
850              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
851 
852              IF ( sec%llstrpond ) THEN
853                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp
854                 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
857              ENDIF
858   
859           ENDDO !end of loop on the level
860
861#if defined key_lim2 || defined key_lim3
862
863           !ICE CASE   
864           !------------
865           IF( sec%ll_ice_section )THEN
866              SELECT CASE (sec%direction(jseg))
867              CASE(0)
868                 zumid_ice = 0
869                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
870              CASE(1)
871                 zumid_ice = 0
872                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
873              CASE(2)
874                 zvmid_ice = 0
875                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
876              CASE(3)
877                 zvmid_ice = 0
878                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
879              END SELECT
880   
881              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
882
883#if defined key_lim2   
884              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   & 
885                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
886                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
887                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
888              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   & 
889                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
890#endif
891#if defined key_lim3
892              DO jl=1,jpl
893                 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*     &
894                                   a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * &
895                                  ( ht_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  &
896                                    ht_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
897                                   
898                 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
899                                   a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
900              ENDDO
901#endif
902   
903           ENDIF !end of ice case
904#endif
905 
906        ENDDO !end of loop on the segment
907
908     ENDIF !end of sec%nb_point =0 case
909     !
910  END SUBROUTINE transport
911 
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 
1126  SUBROUTINE dia_dct_sum(sec,jsec) 
1127     !!-------------------------------------------------------------
1128     !! Purpose: Average the transport over nn_dctwri time steps 
1129     !! and sum over the density/salinity/temperature/depth classes
1130     !!
1131     !! Method:   Sum over relevant grid cells to obtain values 
1132     !!           for each class
1133     !!              There are several loops:                 
1134     !!              loop on the segment between 2 nodes
1135     !!              loop on the level jk
1136     !!              loop on the density/temperature/salinity/level classes
1137     !!              test on the density/temperature/salinity/level
1138     !!
1139     !!  Note:    Transport through a given section is equal to the sum of transports
1140     !!           computed on each proc.
1141     !!           On each proc,transport is equal to the sum of transport computed through
1142     !!           segments linking each point of sec%listPoint  with the next one.   
1143     !!
1144     !!-------------------------------------------------------------
1145     !! * arguments
1146     TYPE(SECTION),INTENT(INOUT) :: sec 
1147     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1148 
1149     TYPE(POINT_SECTION) :: k 
1150     INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes 
1151     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1152     !!-------------------------------------------------------------
1153 
1154     !! Sum the relevant segments to obtain values for each class
1155     IF(sec%nb_point .NE. 0)THEN   
1156 
1157        !--------------------------------------!
1158        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1159        !--------------------------------------!
1160        DO jseg=1,MAX(sec%nb_point-1,0) 
1161           
1162           !-------------------------------------------------------------------------------------------
1163           ! Select the appropriate coordinate for computing the velocity of the segment
1164           !
1165           !                      CASE(0)                                    Case (2)
1166           !                      -------                                    --------
1167           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
1168           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
1169           !                                                                            |
1170           !                                                                            |
1171           !                                                                            |
1172           !                      Case (3)                                            U(i,j)
1173           !                      --------                                              |
1174           !                                                                            |
1175           !  listPoint(jseg+1) F(i,j+1)                                                |
1176           !                        |                                                   |
1177           !                        |                                                   |
1178           !                        |                                 listPoint(jseg+1) F(i,j-1)
1179           !                        |                                             
1180           !                        |                                             
1181           !                     U(i,j+1)                                             
1182           !                        |                                       Case(1)     
1183           !                        |                                       ------       
1184           !                        |                                             
1185           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1186           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
1187           ! listPoint(jseg)     F(i,j)
1188           
1189           !-------------------------------------------------------------------------------------------
1190 
1191           SELECT CASE( sec%direction(jseg) ) 
1192           CASE(0)  ;   k = sec%listPoint(jseg) 
1193           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
1194           CASE(2)  ;   k = sec%listPoint(jseg) 
1195           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
1196           END SELECT 
1197 
1198           !---------------------------|
1199           !     LOOP ON THE LEVEL     |
1200           !---------------------------|
1201           !Sum of the transport on the vertical 
1202           DO jk=1,mbathy(k%I,k%J) 
1203 
1204              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1205              SELECT CASE( sec%direction(jseg) ) 
1206              CASE(0,1) 
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) 
1211
1212              CASE(2,3) 
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 ) 
1217                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1218              END SELECT
1219 
1220              zfsdep= fsdept(k%I,k%J,jk) 
1221 
1222              !-------------------------------
1223              !  LOOP ON THE DENSITY CLASSES |
1224              !-------------------------------
1225              !The computation is made for each density/temperature/salinity/depth class
1226              DO jclass=1,MAX(1,sec%nb_class-1) 
1227 
1228                 !----------------------------------------------!
1229                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 
1230                 !----------------------------------------------!
1231
1232                 IF ( (                                                    & 
1233                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
1234                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
1235                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
1236 
1237                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
1238                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
1239                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
1240 
1241                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
1242                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
1243                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
1244 
1245                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
1246                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
1247                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
1248 
1249                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
1250                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
1251                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
1252                                                                   ))   THEN 
1253 
1254                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
1255                    !----------------------------------------------------------------------------
1256                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN 
1257                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
1258                    ELSE
1259                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
1260                    ENDIF
1261                    IF( sec%llstrpond )THEN
1262 
1263                       IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN
1264                          sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk) 
1265                       ELSE
1266                          sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk) 
1267                       ENDIF
1268 
1269                       IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN
1270                          sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk) 
1271                       ELSE
1272                          sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk) 
1273                       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
1286 
1287                    ELSE
1288                       sec%transport( 3,jclass) = 0._wp 
1289                       sec%transport( 4,jclass) = 0._wp 
1290                       sec%transport( 5,jclass) = 0._wp 
1291                       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
1296                    ENDIF
1297 
1298                 ENDIF ! end of test if point is in class
1299   
1300              ENDDO ! end of loop on the classes
1301 
1302           ENDDO ! loop over jk
1303 
1304#if defined key_lim2 || defined key_lim3 
1305 
1306           !ICE CASE     
1307           IF( sec%ll_ice_section )THEN
1308 
1309              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
1310                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)*1.E-6 
1311              ELSE
1312                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)*1.E-6 
1313              ENDIF
1314 
1315              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
1316                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)*1.E-6 
1317              ELSE
1318                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)*1.E-6 
1319              ENDIF
1320 
1321           ENDIF !end of ice case
1322#endif
1323 
1324        ENDDO !end of loop on the segment
1325 
1326     ELSE  !if sec%nb_point =0
1327        sec%transport(1:2,:)=0. 
1328        IF (sec%llstrpond) sec%transport(3:10,:)=0. 
1329        IF (sec%ll_ice_section) sec%transport(11:14,:)=0. 
1330     ENDIF !end of sec%nb_point =0 case
1331 
1332  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
1544 
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
1721  SUBROUTINE dia_dct_wri(kt,ksec,sec)
1722     !!-------------------------------------------------------------
1723     !! Write transport output in numdct
1724     !!
1725     !! Purpose: Write  transports in ascii files
1726     !!
1727     !! Method:
1728     !!        1. Write volume transports in "volume_transport"
1729     !!           Unit: Sv : area * Velocity / 1.e6
1730     !!
1731     !!        2. Write heat transports in "heat_transport"
1732     !!           Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
1733     !!
1734     !!        3. Write salt transports in "salt_transport"
1735     !!           Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9
1736     !!
1737     !!-------------------------------------------------------------
1738     !!arguments
1739     INTEGER, INTENT(IN)          :: kt          ! time-step
1740     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
1741     INTEGER ,INTENT(IN)          :: ksec        ! section number
1742
1743     !!local declarations
1744     INTEGER               :: jclass             ! Dummy loop
1745     CHARACTER(len=2)      :: classe             ! Classname
1746     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
1747     REAL(wp)              :: zslope             ! section's slope coeff
1748     !
1749     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
1750     !!-------------------------------------------------------------
1751     CALL wrk_alloc(nb_type_class , zsumclasses ) 
1752
1753     zsumclasses(:)=0._wp
1754     zslope = sec%slopeSection       
1755
1756 
1757     DO jclass=1,MAX(1,sec%nb_class-1)
1758
1759        classe   = 'N       '
1760        zbnd1   = 0._wp
1761        zbnd2   = 0._wp
1762        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
1763
1764   
1765        !insitu density classes transports
1766        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
1767            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
1768           classe = 'DI       '
1769           zbnd1 = sec%zsigi(jclass)
1770           zbnd2 = sec%zsigi(jclass+1)
1771        ENDIF
1772        !potential density classes transports
1773        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
1774            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
1775           classe = 'DP      '
1776           zbnd1 = sec%zsigp(jclass)
1777           zbnd2 = sec%zsigp(jclass+1)
1778        ENDIF
1779        !depth classes transports
1780        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
1781            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
1782           classe = 'Z       '
1783           zbnd1 = sec%zlay(jclass)
1784           zbnd2 = sec%zlay(jclass+1)
1785        ENDIF
1786        !salinity classes transports
1787        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
1788            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
1789           classe = 'S       '
1790           zbnd1 = sec%zsal(jclass)
1791           zbnd2 = sec%zsal(jclass+1)   
1792        ENDIF
1793        !temperature classes transports
1794        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
1795            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
1796           classe = 'T       '
1797           zbnd1 = sec%ztem(jclass)
1798           zbnd2 = sec%ztem(jclass+1)
1799        ENDIF
1800                 
1801        !write volume transport per class
1802        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1803                              jclass,classe,zbnd1,zbnd2,&
1804                              sec%transport(1,jclass),sec%transport(2,jclass), &
1805                              sec%transport(1,jclass)+sec%transport(2,jclass)
1806
1807        IF( sec%llstrpond )THEN
1808
1809           !write heat transport per class:
1810           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  &
1811                              jclass,classe,zbnd1,zbnd2,&
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
1814           !write salt transport per class
1815           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  &
1816                              jclass,classe,zbnd1,zbnd2,&
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
1819        ENDIF
1820
1821     ENDDO
1822
1823     zbnd1 = 0._wp
1824     zbnd2 = 0._wp
1825     jclass=0
1826
1827     !write total volume transport
1828     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1829                           jclass,"total",zbnd1,zbnd2,&
1830                           zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
1831
1832     IF( sec%llstrpond )THEN
1833
1834        !write total heat transport
1835        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1836                           jclass,"total",zbnd1,zbnd2,&
1837                           zsumclasses(7)*1.e-15,zsumclasses(8)*1.e-15,&
1838                           (zsumclasses(7)+zsumclasses(8) )*1.e-15
1839        !write total salt transport
1840        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1841                           jclass,"total",zbnd1,zbnd2,&
1842                           zsumclasses(9)*1.e-9,zsumclasses(10)*1.e-9,&
1843                           (zsumclasses(9)+zsumclasses(10))*1.e-9
1844     ENDIF
1845
1846     
1847     IF ( sec%ll_ice_section) THEN
1848        !write total ice volume transport
1849        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1850                              jclass,"ice_vol",zbnd1,zbnd2,&
1851                              sec%transport(11,1),sec%transport(12,1),&
1852                              sec%transport(11,1)+sec%transport(12,1)
1853        !write total ice surface transport
1854        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1855                              jclass,"ice_surf",zbnd1,zbnd2,&
1856                              sec%transport(13,1),sec%transport(14,1), &
1857                              sec%transport(13,1)+sec%transport(14,1) 
1858     ENDIF
1859                                             
1860118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
1861119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
1862
1863     CALL wrk_dealloc(nb_type_class , zsumclasses ) 
1864  END SUBROUTINE dia_dct_wri
1865
1866   PURE  FUNCTION interp(ki, kj, kk, cd_point,var) 
1867  !!----------------------------------------------------------------------
1868  !!
1869  !!   Purpose: compute temperature/salinity/density at U-point or V-point
1870  !!   --------
1871  !!
1872  !!   Method:
1873  !!   ------
1874  !!
1875  !!   ====> full step and partial step
1876  !!
1877  !!
1878  !!    |    I          |    I+1           |    Z=temperature/salinity/density at U-poinT
1879  !!    |               |                  |
1880  !!  ----------------------------------------  1. Veritcal interpolation: compute zbis
1881  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
1882  !!    |               |                  |       zbis =
1883  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
1884  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
1885  !!    |               |                  |
1886  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
1887  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K) 
1888  !!    |     .         |                  |
1889  !!    |     .         |                  |       interp = ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)
1890  !!    |     .         |                  |
1891  !!  ------------------------------------------
1892  !!    |     .         |                  |
1893  !!    |     .         |                  |
1894  !!    |     .         |                  |
1895  !!K   |    zbis.......U...ptab(I+1,J,K)  |
1896  !!    |     .         |                  |
1897  !!    | ptab(I,J,K)   |                  |
1898  !!    |               |------------------|
1899  !!    |               | partials         |
1900  !!    |               |  steps           |
1901  !!  -------------------------------------------
1902  !!    <----zet1------><----zet2--------->
1903  !!
1904  !!
1905  !!   ====>  s-coordinate
1906  !!     
1907  !!    |                |                  |   1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
1908  !!    |                |                  |      Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
1909  !!    |                | ptab(I+1,J,K)    |
1910  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point
1911  !!    |                |      ^           |   
1912  !!    |                |      | zdep2     |   
1913  !!    |                |      |           |   
1914  !!    |       ^        U      v           |
1915  !!    |       |        |                  |
1916  !!    |       | zdep1  |                  |   
1917  !!    |       v        |                  |
1918  !!    |      T1        |                  |
1919  !!    | ptab(I,J,K)    |                  |
1920  !!    |                |                  |
1921  !!    |                |                  |
1922  !!
1923  !!    <----zet1--------><----zet2--------->
1924  !!
1925  !!----------------------------------------------------------------------
1926  !*arguments
1927  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
1928  INTEGER, INTENT(IN)                          :: var   !  which variable
1929  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
1930  REAL(wp)                                     :: interp       ! interpolated variable
1931
1932  !*local declations
1933  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
1934  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real
1935  REAL(wp):: zet1, zet2                                        ! weight for interpolation
1936  REAL(wp):: zdep1,zdep2                                       ! differences of depth
1937  REAL(wp):: zmsk                                              ! mask value
1938  !!----------------------------------------------------------------------
1939
1940 
1941
1942  IF( cd_point=='U' )THEN
1943     ii1 = ki    ; ij1 = kj 
1944     ii2 = ki+1  ; ij2 = kj 
1945
1946     zet1=e1t(ii1,ij1)
1947     zet2=e1t(ii2,ij2)
1948     zmsk=umask(ii1,ij1,kk)
1949 
1950
1951  ELSE ! cd_point=='V'
1952     ii1 = ki    ; ij1 = kj 
1953     ii2 = ki    ; ij2 = kj+1 
1954
1955     zet1=e2t(ii1,ij1)
1956     zet2=e2t(ii2,ij2)
1957     zmsk=vmask(ii1,ij1,kk)
1958
1959  ENDIF
1960
1961  IF( ln_sco )THEN   ! s-coordinate case
1962
1963     zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2 
1964     zdep1 = fsdept(ii1,ij1,kk) - zdepu
1965     zdep2 = fsdept(ii2,ij2,kk) - zdepu
1966
1967     ! weights
1968     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1969     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1970 
1971     ! result
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
1978
1979  ELSE       ! full step or partial step case
1980
1981#if defined key_vvl
1982
1983     ze3t  = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk) 
1984     zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk)
1985     zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk)
1986
1987#else
1988
1989     ze3t  = fse3t(ii2,ij2,kk)   - fse3t(ii1,ij1,kk) 
1990     zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
1991     zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
1992
1993#endif
1994
1995     IF(kk .NE. 1)THEN
1996
1997        IF( ze3t >= 0. )THEN 
1998           ! zbis
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
2013           ! result
2014        ELSE
2015           ! zbis
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
2030        ENDIF   
2031
2032     ELSE
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
2043     ENDIF
2044
2045  ENDIF
2046
2047  END FUNCTION interp
2048
2049#else
2050   !!----------------------------------------------------------------------
2051   !!   Default option :                                       Dummy module
2052   !!----------------------------------------------------------------------
2053   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
2054   PUBLIC 
2055   !! $Id$
2056CONTAINS
2057
2058   SUBROUTINE dia_dct_init          ! Dummy routine
2059      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
2060   END SUBROUTINE dia_dct_init
2061
2062   SUBROUTINE dia_dct( kt )         ! Dummy routine
2063      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
2064      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
2065   END SUBROUTINE dia_dct
2066#endif
2067
2068END MODULE diadct
Note: See TracBrowser for help on using the repository browser.