source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

  • Property svn:executable set to *
File size: 59.9 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
62#include "domzgr_substitute.h90"
63
64  !! * Shared module variables
65  LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
66
67  !! * Module variables
68  INTEGER :: nn_dct        ! Frequency of computation
69  INTEGER :: nn_dctwri     ! Frequency of output
70  INTEGER :: nn_secdebug   ! Number of the section to debug
71   
72  INTEGER, PARAMETER :: nb_class_max  = 10
73  INTEGER, PARAMETER :: nb_sec_max    = 150
74  INTEGER, PARAMETER :: nb_point_max  = 2000
75  INTEGER, PARAMETER :: nb_type_class = 10
76  INTEGER, PARAMETER :: nb_3d_vars    = 3 
77  INTEGER, PARAMETER :: nb_2d_vars    = 2 
78  INTEGER            :: nb_sec 
79
80  TYPE POINT_SECTION
81     INTEGER :: I,J
82  END TYPE POINT_SECTION
83
84  TYPE COORD_SECTION
85     REAL(wp) :: lon,lat
86  END TYPE COORD_SECTION
87
88  TYPE SECTION
89     CHARACTER(len=60)                            :: name              ! name of the sec
90     LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
91                                                                       ! heat transports
92     LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
93     LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
94     TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
95     INTEGER                                      :: nb_class          ! number of boundaries for density classes
96     INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
97     CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
98     REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
99                                                     zsigp           ,&! potential density classes    (99 if you don't want)
100                                                     zsal            ,&! salinity classes   (99 if you don't want)
101                                                     ztem            ,&! temperature classes(99 if you don't want)
102                                                     zlay              ! level classes      (99 if you don't want)
103     REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
104     REAL(wp)                                         :: slopeSection  ! slope of the section
105     INTEGER                                          :: nb_point      ! number of points in the section
106     TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
107  END TYPE SECTION
108
109  TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
110 
111  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
112  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d 
113
114   !! $Id$
115CONTAINS
116
117 
118  INTEGER FUNCTION diadct_alloc() 
119     !!----------------------------------------------------------------------
120     !!                   ***  FUNCTION diadct_alloc  ***
121     !!----------------------------------------------------------------------
122     INTEGER :: ierr(2) 
123     !!----------------------------------------------------------------------
124
125     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) ) 
126     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) ) 
127
128     diadct_alloc = MAXVAL( ierr ) 
129     IF( diadct_alloc /= 0 )   CALL ctl_warn('diadct_alloc: failed to allocate arrays') 
130 
131  END FUNCTION diadct_alloc 
132
133  SUBROUTINE dia_dct_init
134     !!---------------------------------------------------------------------
135     !!               ***  ROUTINE diadct  *** 
136     !!
137     !!  ** Purpose: Read the namelist parameters
138     !!              Open output files
139     !!
140     !!---------------------------------------------------------------------
141     NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
142     INTEGER  ::   ios                 ! Local integer output status for namelist read
143
144     IF( nn_timing == 1 )   CALL timing_start('dia_dct_init')
145
146     REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections
147     READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
148901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp )
149
150     REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
151     READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
152902  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp )
153     IF(lwm) WRITE ( numond, namdct )
154
155     IF( lwp ) THEN
156        WRITE(numout,*) " "
157        WRITE(numout,*) "diadct_init: compute transports through sections "
158        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
159        WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct
160        WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri
161
162        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
163                                            WRITE(numout,*)"       Debug section number: ", nn_secdebug 
164        ELSE IF ( nn_secdebug ==  0 )THEN ; WRITE(numout,*)"       No section to debug"
165        ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)"       Debug all sections"
166        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
167        ENDIF
168
169        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
170          &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
171
172     ENDIF
173
174     !Read section_ijglobal.diadct
175     CALL readsec
176
177     !open output file
178     IF( lwm ) THEN
179        CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
180        CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
181        CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
182     ENDIF
183
184     ! Initialise arrays to zero
185     transports_3d(:,:,:,:)=0.0 
186     transports_2d(:,:,:)  =0.0 
187
188     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init')
189     !
190  END SUBROUTINE dia_dct_init
191 
192 
193  SUBROUTINE dia_dct(kt)
194     !!---------------------------------------------------------------------
195     !!               ***  ROUTINE diadct  *** 
196     !!
197     !!  Purpose :: Compute section transports and write it in numdct files
198     !!   
199     !!  Method  :: All arrays initialised to zero in dct_init
200     !!             Each nn_dct time step call subroutine 'transports' for
201     !!               each section to sum the transports over each grid cell.
202     !!             Each nn_dctwri time step:
203     !!               Divide the arrays by the number of summations to gain
204     !!               an average value
205     !!               Call dia_dct_sum to sum relevant grid boxes to obtain
206     !!               totals for each class (density, depth, temp or sal)
207     !!               Call dia_dct_wri to write the transports into file
208     !!               Reinitialise all relevant arrays to zero
209     !!---------------------------------------------------------------------
210     !! * Arguments
211     INTEGER,INTENT(IN)        ::kt
212
213     !! * Local variables
214     INTEGER             :: jsec,            &! loop on sections
215                            itotal            ! nb_sec_max*nb_type_class*nb_class_max
216     LOGICAL             :: lldebug =.FALSE.  ! debug a section 
217     
218     INTEGER , DIMENSION(1)             :: ish   ! tmp array for mpp_sum
219     INTEGER , DIMENSION(3)             :: ish2  !   "
220     REAL(wp), POINTER, DIMENSION(:)    :: zwork !   " 
221     REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   "
222
223     !!---------------------------------------------------------------------   
224     IF( nn_timing == 1 )   CALL timing_start('dia_dct')
225
226     IF( lk_mpp )THEN
227        itotal = nb_sec_max*nb_type_class*nb_class_max
228        CALL wrk_alloc( itotal                                , zwork ) 
229        CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum  )
230     ENDIF   
231 
232     ! Initialise arrays
233     zwork(:) = 0.0 
234     zsum(:,:,:) = 0.0
235
236     IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
237         WRITE(numout,*) " "
238         WRITE(numout,*) "diadct: compute transport"
239         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
240         WRITE(numout,*) "nb_sec = ",nb_sec
241     ENDIF
242
243 
244     ! Compute transport and write only at nn_dctwri
245     IF( MOD(kt,nn_dct)==0 ) THEN
246
247        DO jsec=1,nb_sec
248
249           !debug this section computing ?
250           lldebug=.FALSE.
251           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. 
252
253           !Compute transport through section 
254           CALL transport(secs(jsec),lldebug,jsec) 
255
256        ENDDO
257             
258        IF( MOD(kt,nn_dctwri)==0 )THEN
259
260           IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: average transports and write at kt = ",kt         
261 
262           !! divide arrays by nn_dctwri/nn_dct to obtain average
263           transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 
264           transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct) 
265 
266           ! Sum over each class
267           DO jsec=1,nb_sec 
268              CALL dia_dct_sum(secs(jsec),jsec) 
269           ENDDO 
270
271           !Sum on all procs
272           IF( lk_mpp )THEN
273              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
274              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
275              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
276              zwork(:)= RESHAPE(zsum(:,:,:), ish )
277              CALL mpp_sum(zwork, ish(1))
278              zsum(:,:,:)= RESHAPE(zwork,ish2)
279              DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
280           ENDIF
281
282           !Write the transport
283           DO jsec=1,nb_sec
284
285              IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
286           
287              !nullify transports values after writing
288              transports_3d(:,jsec,:,:)=0.
289              transports_2d(:,jsec,:  )=0.
290              secs(jsec)%transport(:,:)=0. 
291
292           ENDDO
293
294        ENDIF
295
296     ENDIF
297
298     IF( lk_mpp )THEN
299        itotal = nb_sec_max*nb_type_class*nb_class_max
300        CALL wrk_dealloc( itotal                                , zwork ) 
301        CALL wrk_dealloc( nb_sec_max,nb_type_class,nb_class_max , zsum  )
302     ENDIF   
303
304     IF( nn_timing == 1 )   CALL timing_stop('dia_dct')
305     !
306  END SUBROUTINE dia_dct
307
308  SUBROUTINE readsec 
309     !!---------------------------------------------------------------------
310     !!               ***  ROUTINE readsec  ***
311     !!
312     !!  ** Purpose:
313     !!            Read a binary file(section_ijglobal.diadct)
314     !!            generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
315     !!
316     !!
317     !!---------------------------------------------------------------------
318     !! * Local variables
319     INTEGER :: iptglo , iptloc                               ! Global and local number of points for a section
320     INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2  ! temporary  integer
321     INTEGER :: jsec, jpt                                     ! dummy loop indices
322
323     INTEGER, DIMENSION(2) :: icoord 
324     CHARACTER(len=160)    :: clname                          !filename
325     CHARACTER(len=200)    :: cltmp
326     CHARACTER(len=200)    :: clformat                        !automatic format
327     TYPE(POINT_SECTION),DIMENSION(nb_point_max)  ::coordtemp !contains listpoints coordinates
328                                                              !read in the file
329     INTEGER, POINTER, DIMENSION(:) :: directemp              !contains listpoints directions
330                                                              !read in the files
331     LOGICAL :: llbon                                       ,&!local logical
332                lldebug                                       !debug the section
333     !!-------------------------------------------------------------------------------------
334     CALL wrk_alloc( nb_point_max, directemp )
335
336     !open input file
337     !---------------
338     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
339 
340     !---------------
341     !Read input file
342     !---------------
343     
344     DO jsec=1,nb_sec_max      !loop on the nb_sec sections
345
346        IF ( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) ) &
347           & WRITE(numout,*)'debuging for section number: ',jsec 
348
349        !initialization
350        !---------------
351        secs(jsec)%name=''
352        secs(jsec)%llstrpond    = .FALSE.  ; secs(jsec)%ll_ice_section = .FALSE.
353        secs(jsec)%ll_date_line = .FALSE.  ; secs(jsec)%nb_class       = 0
354        secs(jsec)%zsigi        = 99._wp   ; secs(jsec)%zsigp          = 99._wp
355        secs(jsec)%zsal         = 99._wp   ; secs(jsec)%ztem           = 99._wp
356        secs(jsec)%zlay         = 99._wp         
357        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0
358
359        !read section's number / name / computing choices / classes / slopeSection / points number
360        !-----------------------------------------------------------------------------------------
361        READ(numdct_in,iostat=iost)isec
362        IF (iost .NE. 0 )EXIT !end of file
363        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
364        IF( jsec .NE. isec )  CALL ctl_stop( cltmp )
365
366        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec 
367
368        READ(numdct_in)secs(jsec)%name
369        READ(numdct_in)secs(jsec)%llstrpond
370        READ(numdct_in)secs(jsec)%ll_ice_section
371        READ(numdct_in)secs(jsec)%ll_date_line
372        READ(numdct_in)secs(jsec)%coordSec
373        READ(numdct_in)secs(jsec)%nb_class
374        READ(numdct_in)secs(jsec)%zsigi
375        READ(numdct_in)secs(jsec)%zsigp
376        READ(numdct_in)secs(jsec)%zsal
377        READ(numdct_in)secs(jsec)%ztem
378        READ(numdct_in)secs(jsec)%zlay
379        READ(numdct_in)secs(jsec)%slopeSection
380        READ(numdct_in)iptglo
381
382        !debug
383        !-----
384
385        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
386         
387            WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))' 
388
389            WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name)
390            WRITE(numout,*)       "      Compute heat and salt transport ? ",secs(jsec)%llstrpond
391            WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section
392            WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line
393            WRITE(numout,*)       "      Slope section :                   ",secs(jsec)%slopeSection
394            WRITE(numout,*)       "      Number of points in the section:  ",iptglo
395            WRITE(numout,*)       "      Number of classes                 ",secs(jsec)%nb_class
396            WRITE(numout,clformat)"      Insitu density classes :          ",secs(jsec)%zsigi
397            WRITE(numout,clformat)"      Potential density classes :       ",secs(jsec)%zsigp
398            WRITE(numout,clformat)"      Salinity classes :                ",secs(jsec)%zsal
399            WRITE(numout,clformat)"      Temperature classes :             ",secs(jsec)%ztem
400            WRITE(numout,clformat)"      Depth classes :                   ",secs(jsec)%zlay
401        ENDIF               
402
403        IF( iptglo .NE. 0 )THEN
404             
405           !read points'coordinates and directions
406           !--------------------------------------
407           coordtemp(:) = POINT_SECTION(0,0) !list of points read
408           directemp(:) = 0                  !value of directions of each points
409           DO jpt=1,iptglo
410              READ(numdct_in)i1,i2
411              coordtemp(jpt)%I = i1 
412              coordtemp(jpt)%J = i2
413           ENDDO
414           READ(numdct_in)directemp(1:iptglo)
415   
416           !debug
417           !-----
418           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
419              WRITE(numout,*)"      List of points in global domain:"
420              DO jpt=1,iptglo
421                 WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt),directemp(jpt)
422              ENDDO                 
423           ENDIF
424 
425           !Now each proc selects only points that are in its domain:
426           !--------------------------------------------------------
427           iptloc = 0                    !initialize number of points selected
428           DO jpt=1,iptglo               !loop on listpoint read in the file
429                   
430              iiglo=coordtemp(jpt)%I          ! global coordinates of the point
431              ijglo=coordtemp(jpt)%J          !  "
432
433              IF( iiglo==jpidta .AND. nimpp==1 ) iiglo = 2
434
435              iiloc=iiglo-jpizoom+1-nimpp+1   ! local coordinates of the point
436              ijloc=ijglo-jpjzoom+1-njmpp+1   !  "
437
438              !verify if the point is on the local domain:(1,nlei)*(1,nlej)
439              IF( iiloc .GE. 1 .AND. iiloc .LE. nlei .AND. &
440                  ijloc .GE. 1 .AND. ijloc .LE. nlej       )THEN
441                 iptloc = iptloc + 1                                                 ! count local points
442                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
443                 secs(jsec)%direction(iptloc) = directemp(jpt)                       ! store local direction
444              ENDIF
445
446           ENDDO
447     
448           secs(jsec)%nb_point=iptloc !store number of section's points
449
450           !debug
451           !-----
452           IF(   lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
453              WRITE(numout,*)"      List of points selected by the proc:"
454              DO jpt = 1,iptloc
455                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
456                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
457                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
458              ENDDO
459           ENDIF
460
461              IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
462              DO jpt = 1,iptloc
463                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
464                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
465              ENDDO
466              ENDIF
467
468           !remove redundant points between processors
469           !------------------------------------------
470           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
471           IF( iptloc .NE. 0 )THEN
472              CALL removepoints(secs(jsec),'I','top_list',lldebug)
473              CALL removepoints(secs(jsec),'I','bot_list',lldebug)
474              CALL removepoints(secs(jsec),'J','top_list',lldebug)
475              CALL removepoints(secs(jsec),'J','bot_list',lldebug)
476           ENDIF
477           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
478              DO jpt = 1,secs(jsec)%nb_point
479                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
480                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
481              ENDDO
482           ENDIF
483
484           !debug
485           !-----
486           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
487              WRITE(numout,*)"      List of points after removepoints:"
488              iptloc = secs(jsec)%nb_point
489              DO jpt = 1,iptloc
490                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
491                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
492                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
493                 CALL FLUSH(numout)
494              ENDDO
495           ENDIF
496
497        ELSE  ! iptglo = 0
498           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )&
499              WRITE(numout,*)'   No points for this section.'
500        ENDIF
501
502     ENDDO !end of the loop on jsec
503 
504     nb_sec = jsec-1   !number of section read in the file
505
506     CALL wrk_dealloc( nb_point_max, directemp )
507     !
508  END SUBROUTINE readsec
509
510  SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
511     !!---------------------------------------------------------------------------
512     !!             *** function removepoints
513     !!
514     !!   ** Purpose :: Remove points which are common to 2 procs
515     !!
516     !----------------------------------------------------------------------------
517     !! * arguments
518     TYPE(SECTION),INTENT(INOUT) :: sec
519     CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
520     CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
521     LOGICAL,INTENT(IN)          :: ld_debug                     
522
523     !! * Local variables
524     INTEGER :: iextr         ,& !extremity of listpoint that we verify
525                iind          ,& !coord     of listpoint that we verify
526                itest         ,& !indice value of the side of the domain
527                                 !where points could be redundant
528                isgn          ,& ! isgn= 1 : scan listpoint from start to end
529                                 ! isgn=-1 : scan listpoint from end to start
530                istart,iend      !first and last points selected in listpoint
531     INTEGER :: jpoint           !loop on list points
532     INTEGER, POINTER, DIMENSION(:)   :: idirec !contains temporary sec%direction
533     INTEGER, POINTER, DIMENSION(:,:) :: icoord !contains temporary sec%listpoint
534     !----------------------------------------------------------------------------
535     CALL wrk_alloc(    nb_point_max, idirec )
536     CALL wrk_alloc( 2, nb_point_max, icoord )
537
538     IF( ld_debug )WRITE(numout,*)'      -------------------------'
539     IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
540
541     !iextr=extremity of list_point that we verify
542     IF      ( cdextr=='bot_list' )THEN ; iextr=1            ; isgn=1
543     ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
544     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdextr")
545     ENDIF
546 
547     !which coordinate shall we verify ?
548     IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1
549     ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2
550     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind") 
551     ENDIF
552
553     IF( ld_debug )THEN
554        WRITE(numout,*)'      case: coord/list extr/domain side'
555        WRITE(numout,*)'      ', cdind,' ',cdextr,' ',itest
556        WRITE(numout,*)'      Actual number of points: ',sec%nb_point
557     ENDIF
558
559     icoord(1,1:nb_point_max) = sec%listPoint%I
560     icoord(2,1:nb_point_max) = sec%listPoint%J
561     idirec                   = sec%direction
562     sec%listPoint            = POINT_SECTION(0,0)
563     sec%direction            = 0
564
565     jpoint=iextr+isgn
566     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
567         IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
568         ELSE                                                                               ; EXIT
569         ENDIF
570     ENDDO 
571
572     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
573     ELSE                        ; istart=1        ; iend=jpoint+1
574     ENDIF
575
576     sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
577     sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
578     sec%direction(1:1+iend-istart)   = idirec(istart:iend)
579     sec%nb_point                     = iend-istart+1
580     
581     IF( ld_debug )THEN
582        WRITE(numout,*)'      Number of points after removepoints :',sec%nb_point
583        WRITE(numout,*)'      sec%direction after removepoints :',sec%direction(1:sec%nb_point)
584     ENDIF
585
586     CALL wrk_dealloc(    nb_point_max, idirec )
587     CALL wrk_dealloc( 2, nb_point_max, icoord )
588  END SUBROUTINE removepoints
589
590  SUBROUTINE transport(sec,ld_debug,jsec)
591     !!-------------------------------------------------------------------------------------------
592     !!                     ***  ROUTINE transport  ***
593     !!
594     !!  Purpose ::  Compute the transport for each point in a section
595     !!
596     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
597     !!              Be aware :           
598     !!              One section is a sum of segments
599     !!              One segment is defined by 2 consecutive points in sec%listPoint
600     !!              All points of sec%listPoint are positioned on the F-point of the cell
601     !!
602     !!              There are two loops:                 
603     !!              loop on the segment between 2 nodes
604     !!              loop on the level jk !!
605     !!
606     !!  Output  ::  Arrays containing the volume,density,heat,salt transports for each i
607     !!              point in a section, summed over each nn_dct.
608     !!
609     !!-------------------------------------------------------------------------------------------
610     !! * Arguments
611     TYPE(SECTION),INTENT(INOUT) :: sec
612     LOGICAL      ,INTENT(IN)    :: ld_debug
613     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
614   
615     !! * Local variables
616     INTEGER             :: jk, jseg, jclass,jl,                 &!loop on level/segment/classes/ice categories
617                            isgnu, isgnv                          !
618     REAL(wp)            :: zumid, zvmid,                        &!U/V velocity on a cell segment
619                            zumid_ice, zvmid_ice,                &!U/V ice velocity
620                            zTnorm                                !transport of velocity through one cell's sides
621     REAL(wp)            :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point
622
623     TYPE(POINT_SECTION) :: k
624     !!--------------------------------------------------------
625
626     IF( ld_debug )WRITE(numout,*)'      Compute transport'
627
628     !---------------------------!
629     !  COMPUTE TRANSPORT        !
630     !---------------------------!
631     IF(sec%nb_point .NE. 0)THEN   
632
633        !----------------------------------------------------------------------------------------------------
634        !Compute sign for velocities:
635        !
636        !convention:
637        !   non horizontal section: direction + is toward left hand of section
638        !       horizontal section: direction + is toward north of section
639        !
640        !
641        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
642        !       ----------------      -----------------     ---------------             --------------
643        !
644        !   isgnv=1         direction +     
645        !  ______         _____             ______                                                   
646        !        |           //|            |                  |                         direction +   
647        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
648        !        |_______  //         ______|    \\            | ---\                        |
649        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
650        !               |             |          __\\|         |                   
651        !               |             |     direction +        |                      isgnv=1                                 
652        !                                                     
653        !----------------------------------------------------------------------------------------------------
654        isgnu = 1
655        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
656        ELSE                                ; isgnv =  1
657        ENDIF
658        IF( sec%slopeSection .GE. 9999. )     isgnv =  1
659
660        IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
661
662        !--------------------------------------!
663        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
664        !--------------------------------------!
665        DO jseg=1,MAX(sec%nb_point-1,0)
666             
667           !-------------------------------------------------------------------------------------------
668           ! Select the appropriate coordinate for computing the velocity of the segment
669           !
670           !                      CASE(0)                                    Case (2)
671           !                      -------                                    --------
672           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
673           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
674           !                                                                            |
675           !                                                                            |
676           !                                                                            |
677           !                      Case (3)                                            U(i,j)
678           !                      --------                                              |
679           !                                                                            |
680           !  listPoint(jseg+1) F(i,j+1)                                                |
681           !                        |                                                   |
682           !                        |                                                   |
683           !                        |                                 listPoint(jseg+1) F(i,j-1)
684           !                        |                                           
685           !                        |                                           
686           !                     U(i,j+1)                                           
687           !                        |                                       Case(1)     
688           !                        |                                       ------     
689           !                        |                                           
690           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
691           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
692           ! listPoint(jseg)     F(i,j)
693           !
694           !-------------------------------------------------------------------------------------------
695
696           SELECT CASE( sec%direction(jseg) )
697           CASE(0)  ;   k = sec%listPoint(jseg)
698           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
699           CASE(2)  ;   k = sec%listPoint(jseg)
700           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
701           END SELECT
702
703           !---------------------------|
704           !     LOOP ON THE LEVEL     |
705           !---------------------------|
706           !Sum of the transport on the vertical 
707           DO jk=1,mbathy(k%I,k%J) 
708 
709              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
710              SELECT CASE( sec%direction(jseg) ) 
711              CASE(0,1) 
712                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
713                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
714                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
715                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
716                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
717              CASE(2,3) 
718                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
719                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
720                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
721                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
722                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
723              END SELECT
724 
725              zfsdep= fsdept(k%I,k%J,jk) 
726 
727              !compute velocity with the correct direction
728              SELECT CASE( sec%direction(jseg) ) 
729              CASE(0,1)   
730                 zumid=0. 
731                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
732              CASE(2,3) 
733                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
734                 zvmid=0. 
735              END SELECT 
736 
737              !zTnorm=transport through one cell;
738              !velocity* cell's length * cell's thickness
739              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
740                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
741
742#if ! defined key_vvl
743              !add transport due to free surface
744              IF( jk==1 )THEN
745                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
746                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
747              ENDIF 
748#endif
749              !COMPUTE TRANSPORT 
750 
751              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
752 
753              IF ( sec%llstrpond ) THEN
754                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp
755                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001
756              ENDIF
757   
758           ENDDO !end of loop on the level
759
760#if defined key_lim2 || defined key_lim3
761
762           !ICE CASE   
763           !------------
764           IF( sec%ll_ice_section )THEN
765              SELECT CASE (sec%direction(jseg))
766              CASE(0)
767                 zumid_ice = 0
768                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
769              CASE(1)
770                 zumid_ice = 0
771                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
772              CASE(2)
773                 zvmid_ice = 0
774                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
775              CASE(3)
776                 zvmid_ice = 0
777                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
778              END SELECT
779   
780              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
781
782#if defined key_lim2   
783              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   & 
784                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
785                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
786                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
787              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   & 
788                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
789#endif
790#if defined key_lim3
791              DO jl=1,jpl
792                 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*     &
793                                   a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * &
794                                  ( ht_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  &
795                                    ht_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
796                                   
797                 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
798                                   a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
799              ENDDO
800#endif
801   
802           ENDIF !end of ice case
803#endif
804 
805        ENDDO !end of loop on the segment
806
807     ENDIF !end of sec%nb_point =0 case
808     !
809  END SUBROUTINE transport
810 
811  SUBROUTINE dia_dct_sum(sec,jsec) 
812     !!-------------------------------------------------------------
813     !! Purpose: Average the transport over nn_dctwri time steps 
814     !! and sum over the density/salinity/temperature/depth classes
815     !!
816     !! Method:   Sum over relevant grid cells to obtain values 
817     !!           for each class
818     !!              There are several loops:                 
819     !!              loop on the segment between 2 nodes
820     !!              loop on the level jk
821     !!              loop on the density/temperature/salinity/level classes
822     !!              test on the density/temperature/salinity/level
823     !!
824     !!  Note:    Transport through a given section is equal to the sum of transports
825     !!           computed on each proc.
826     !!           On each proc,transport is equal to the sum of transport computed through
827     !!           segments linking each point of sec%listPoint  with the next one.   
828     !!
829     !!-------------------------------------------------------------
830     !! * arguments
831     TYPE(SECTION),INTENT(INOUT) :: sec 
832     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
833 
834     TYPE(POINT_SECTION) :: k 
835     INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes 
836     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
837     !!-------------------------------------------------------------
838 
839     !! Sum the relevant segments to obtain values for each class
840     IF(sec%nb_point .NE. 0)THEN   
841 
842        !--------------------------------------!
843        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
844        !--------------------------------------!
845        DO jseg=1,MAX(sec%nb_point-1,0) 
846           
847           !-------------------------------------------------------------------------------------------
848           ! Select the appropriate coordinate for computing the velocity of the segment
849           !
850           !                      CASE(0)                                    Case (2)
851           !                      -------                                    --------
852           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
853           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
854           !                                                                            |
855           !                                                                            |
856           !                                                                            |
857           !                      Case (3)                                            U(i,j)
858           !                      --------                                              |
859           !                                                                            |
860           !  listPoint(jseg+1) F(i,j+1)                                                |
861           !                        |                                                   |
862           !                        |                                                   |
863           !                        |                                 listPoint(jseg+1) F(i,j-1)
864           !                        |                                             
865           !                        |                                             
866           !                     U(i,j+1)                                             
867           !                        |                                       Case(1)     
868           !                        |                                       ------       
869           !                        |                                             
870           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
871           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
872           ! listPoint(jseg)     F(i,j)
873           
874           !-------------------------------------------------------------------------------------------
875 
876           SELECT CASE( sec%direction(jseg) ) 
877           CASE(0)  ;   k = sec%listPoint(jseg) 
878           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
879           CASE(2)  ;   k = sec%listPoint(jseg) 
880           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
881           END SELECT 
882 
883           !---------------------------|
884           !     LOOP ON THE LEVEL     |
885           !---------------------------|
886           !Sum of the transport on the vertical 
887           DO jk=1,mbathy(k%I,k%J) 
888 
889              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
890              SELECT CASE( sec%direction(jseg) ) 
891              CASE(0,1) 
892                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
893                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
894                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
895                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
896
897              CASE(2,3) 
898                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
899                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
900                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
901                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
902                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
903              END SELECT
904 
905              zfsdep= fsdept(k%I,k%J,jk) 
906 
907              !-------------------------------
908              !  LOOP ON THE DENSITY CLASSES |
909              !-------------------------------
910              !The computation is made for each density/temperature/salinity/depth class
911              DO jclass=1,MAX(1,sec%nb_class-1) 
912 
913                 !----------------------------------------------!
914                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 
915                 !----------------------------------------------!
916
917                 IF ( (                                                    & 
918                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
919                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
920                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
921 
922                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
923                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
924                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
925 
926                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
927                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
928                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
929 
930                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
931                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
932                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
933 
934                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
935                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
936                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
937                                                                   ))   THEN 
938 
939                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
940                    !----------------------------------------------------------------------------
941                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN 
942                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
943                    ELSE
944                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
945                    ENDIF
946                    IF( sec%llstrpond )THEN
947 
948                       IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN
949                          sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk) 
950                       ELSE
951                          sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk) 
952                       ENDIF
953 
954                       IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN
955                          sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk) 
956                       ELSE
957                          sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk) 
958                       ENDIF
959 
960                    ELSE
961                       sec%transport( 3,jclass) = 0._wp 
962                       sec%transport( 4,jclass) = 0._wp 
963                       sec%transport( 5,jclass) = 0._wp 
964                       sec%transport( 6,jclass) = 0._wp 
965                    ENDIF
966 
967                 ENDIF ! end of test if point is in class
968   
969              ENDDO ! end of loop on the classes
970 
971           ENDDO ! loop over jk
972 
973#if defined key_lim2 || defined key_lim3 
974 
975           !ICE CASE     
976           IF( sec%ll_ice_section )THEN
977 
978              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
979                 sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6 
980              ELSE
981                 sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6 
982              ENDIF
983 
984              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
985                 sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6 
986              ELSE
987                 sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6 
988              ENDIF
989 
990           ENDIF !end of ice case
991#endif
992 
993        ENDDO !end of loop on the segment
994 
995     ELSE  !if sec%nb_point =0
996        sec%transport(1:2,:)=0. 
997        IF (sec%llstrpond) sec%transport(3:6,:)=0. 
998        IF (sec%ll_ice_section) sec%transport(7:10,:)=0. 
999     ENDIF !end of sec%nb_point =0 case
1000 
1001  END SUBROUTINE dia_dct_sum 
1002 
1003  SUBROUTINE dia_dct_wri(kt,ksec,sec)
1004     !!-------------------------------------------------------------
1005     !! Write transport output in numdct
1006     !!
1007     !! Purpose: Write  transports in ascii files
1008     !!
1009     !! Method:
1010     !!        1. Write volume transports in "volume_transport"
1011     !!           Unit: Sv : area * Velocity / 1.e6
1012     !!
1013     !!        2. Write heat transports in "heat_transport"
1014     !!           Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
1015     !!
1016     !!        3. Write salt transports in "salt_transport"
1017     !!           Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9
1018     !!
1019     !!-------------------------------------------------------------
1020     !!arguments
1021     INTEGER, INTENT(IN)          :: kt          ! time-step
1022     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
1023     INTEGER ,INTENT(IN)          :: ksec        ! section number
1024
1025     !!local declarations
1026     INTEGER               :: jclass             ! Dummy loop
1027     CHARACTER(len=2)      :: classe             ! Classname
1028     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
1029     REAL(wp)              :: zslope             ! section's slope coeff
1030     !
1031     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
1032     !!-------------------------------------------------------------
1033     CALL wrk_alloc(nb_type_class , zsumclasses ) 
1034
1035     zsumclasses(:)=0._wp
1036     zslope = sec%slopeSection       
1037
1038 
1039     DO jclass=1,MAX(1,sec%nb_class-1)
1040
1041        classe   = 'N       '
1042        zbnd1   = 0._wp
1043        zbnd2   = 0._wp
1044        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
1045
1046   
1047        !insitu density classes transports
1048        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
1049            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
1050           classe = 'DI       '
1051           zbnd1 = sec%zsigi(jclass)
1052           zbnd2 = sec%zsigi(jclass+1)
1053        ENDIF
1054        !potential density classes transports
1055        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
1056            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
1057           classe = 'DP      '
1058           zbnd1 = sec%zsigp(jclass)
1059           zbnd2 = sec%zsigp(jclass+1)
1060        ENDIF
1061        !depth classes transports
1062        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
1063            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
1064           classe = 'Z       '
1065           zbnd1 = sec%zlay(jclass)
1066           zbnd2 = sec%zlay(jclass+1)
1067        ENDIF
1068        !salinity classes transports
1069        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
1070            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
1071           classe = 'S       '
1072           zbnd1 = sec%zsal(jclass)
1073           zbnd2 = sec%zsal(jclass+1)   
1074        ENDIF
1075        !temperature classes transports
1076        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
1077            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
1078           classe = 'T       '
1079           zbnd1 = sec%ztem(jclass)
1080           zbnd2 = sec%ztem(jclass+1)
1081        ENDIF
1082                 
1083        !write volume transport per class
1084        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1085                              jclass,classe,zbnd1,zbnd2,&
1086                              sec%transport(1,jclass),sec%transport(2,jclass), &
1087                              sec%transport(1,jclass)+sec%transport(2,jclass)
1088
1089        IF( sec%llstrpond )THEN
1090
1091           !write heat transport per class:
1092           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  &
1093                              jclass,classe,zbnd1,zbnd2,&
1094                              sec%transport(3,jclass)*1.e-15,sec%transport(4,jclass)*1.e-15, &
1095                              ( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15
1096           !write salt transport per class
1097           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  &
1098                              jclass,classe,zbnd1,zbnd2,&
1099                              sec%transport(5,jclass)*1.e-9,sec%transport(6,jclass)*1.e-9,&
1100                              (sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9
1101        ENDIF
1102
1103     ENDDO
1104
1105     zbnd1 = 0._wp
1106     zbnd2 = 0._wp
1107     jclass=0
1108
1109     !write total volume transport
1110     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1111                           jclass,"total",zbnd1,zbnd2,&
1112                           zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
1113
1114     IF( sec%llstrpond )THEN
1115
1116        !write total heat transport
1117        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1118                           jclass,"total",zbnd1,zbnd2,&
1119                           zsumclasses(3)*1.e-15,zsumclasses(4)*1.e-15,&
1120                           (zsumclasses(3)+zsumclasses(4) )*1.e-15
1121        !write total salt transport
1122        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1123                           jclass,"total",zbnd1,zbnd2,&
1124                           zsumclasses(5)*1.e-9,zsumclasses(6)*1.e-9,&
1125                           (zsumclasses(5)+zsumclasses(6))*1.e-9
1126     ENDIF
1127
1128     
1129     IF ( sec%ll_ice_section) THEN
1130        !write total ice volume transport
1131        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1132                              jclass,"ice_vol",zbnd1,zbnd2,&
1133                              sec%transport(7,1),sec%transport(8,1),&
1134                              sec%transport(7,1)+sec%transport(8,1)
1135        !write total ice surface transport
1136        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1137                              jclass,"ice_surf",zbnd1,zbnd2,&
1138                              sec%transport(9,1),sec%transport(10,1), &
1139                              sec%transport(9,1)+sec%transport(10,1) 
1140     ENDIF
1141                                             
1142118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
1143119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
1144
1145     CALL wrk_dealloc(nb_type_class , zsumclasses ) 
1146  END SUBROUTINE dia_dct_wri
1147
1148  FUNCTION interp(ki, kj, kk, cd_point, ptab)
1149  !!----------------------------------------------------------------------
1150  !!
1151  !!   Purpose: compute temperature/salinity/density at U-point or V-point
1152  !!   --------
1153  !!
1154  !!   Method:
1155  !!   ------
1156  !!
1157  !!   ====> full step and partial step
1158  !!
1159  !!
1160  !!    |    I          |    I+1           |    Z=temperature/salinity/density at U-poinT
1161  !!    |               |                  |
1162  !!  ----------------------------------------  1. Veritcal interpolation: compute zbis
1163  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
1164  !!    |               |                  |       zbis =
1165  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
1166  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
1167  !!    |               |                  |
1168  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
1169  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K) 
1170  !!    |     .         |                  |
1171  !!    |     .         |                  |       interp = ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)
1172  !!    |     .         |                  |
1173  !!  ------------------------------------------
1174  !!    |     .         |                  |
1175  !!    |     .         |                  |
1176  !!    |     .         |                  |
1177  !!K   |    zbis.......U...ptab(I+1,J,K)  |
1178  !!    |     .         |                  |
1179  !!    | ptab(I,J,K)   |                  |
1180  !!    |               |------------------|
1181  !!    |               | partials         |
1182  !!    |               |  steps           |
1183  !!  -------------------------------------------
1184  !!    <----zet1------><----zet2--------->
1185  !!
1186  !!
1187  !!   ====>  s-coordinate
1188  !!     
1189  !!    |                |                  |   1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
1190  !!    |                |                  |      Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
1191  !!    |                | ptab(I+1,J,K)    |
1192  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point
1193  !!    |                |      ^           |   
1194  !!    |                |      | zdep2     |   
1195  !!    |                |      |           |   
1196  !!    |       ^        U      v           |
1197  !!    |       |        |                  |
1198  !!    |       | zdep1  |                  |   
1199  !!    |       v        |                  |
1200  !!    |      T1        |                  |
1201  !!    | ptab(I,J,K)    |                  |
1202  !!    |                |                  |
1203  !!    |                |                  |
1204  !!
1205  !!    <----zet1--------><----zet2--------->
1206  !!
1207  !!----------------------------------------------------------------------
1208  !*arguments
1209  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
1210  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
1211  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk )
1212  REAL(wp)                                     :: interp       ! interpolated variable
1213
1214  !*local declations
1215  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
1216  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real
1217  REAL(wp):: zet1, zet2                                        ! weight for interpolation
1218  REAL(wp):: zdep1,zdep2                                       ! differences of depth
1219  REAL(wp):: zmsk                                              ! mask value
1220  !!----------------------------------------------------------------------
1221
1222  IF( cd_point=='U' )THEN
1223     ii1 = ki    ; ij1 = kj 
1224     ii2 = ki+1  ; ij2 = kj 
1225
1226     zet1=e1t(ii1,ij1)
1227     zet2=e1t(ii2,ij2)
1228     zmsk=umask(ii1,ij1,kk)
1229 
1230
1231  ELSE ! cd_point=='V'
1232     ii1 = ki    ; ij1 = kj 
1233     ii2 = ki    ; ij2 = kj+1 
1234
1235     zet1=e2t(ii1,ij1)
1236     zet2=e2t(ii2,ij2)
1237     zmsk=vmask(ii1,ij1,kk)
1238
1239  ENDIF
1240
1241  IF( ln_sco )THEN   ! s-coordinate case
1242
1243     zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2 
1244     zdep1 = fsdept(ii1,ij1,kk) - zdepu
1245     zdep2 = fsdept(ii2,ij2,kk) - zdepu
1246
1247     ! weights
1248     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1249     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1250 
1251     ! result
1252     interp = zmsk * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )   
1253
1254
1255  ELSE       ! full step or partial step case
1256
1257#if defined key_vvl
1258
1259     ze3t  = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk) 
1260     zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk)
1261     zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk)
1262
1263#else
1264
1265     ze3t  = fse3t(ii2,ij2,kk)   - fse3t(ii1,ij1,kk) 
1266     zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
1267     zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
1268
1269#endif
1270
1271     IF(kk .NE. 1)THEN
1272
1273        IF( ze3t >= 0. )THEN 
1274           ! zbis
1275           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
1276           ! result
1277            interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
1278        ELSE
1279           ! zbis
1280           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1281           ! result
1282           interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1283        ENDIF   
1284
1285     ELSE
1286        interp = zmsk * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1287     ENDIF
1288
1289  ENDIF
1290
1291
1292  END FUNCTION interp
1293
1294#else
1295   !!----------------------------------------------------------------------
1296   !!   Default option :                                       Dummy module
1297   !!----------------------------------------------------------------------
1298   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
1299   PUBLIC 
1300   !! $Id$
1301CONTAINS
1302
1303   SUBROUTINE dia_dct_init          ! Dummy routine
1304      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
1305   END SUBROUTINE dia_dct_init
1306
1307   SUBROUTINE dia_dct( kt )         ! Dummy routine
1308      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
1309      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
1310   END SUBROUTINE dia_dct
1311#endif
1312
1313END MODULE diadct
Note: See TracBrowser for help on using the repository browser.