New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diadct.F90 in NEMO/branches/UKMO/r8395_restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: NEMO/branches/UKMO/r8395_restart_datestamp/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 10758

Last change on this file since 10758 was 10758, checked in by jcastill, 5 years ago

Remove svn keywords

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