source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 4147

Last change on this file since 4147 was 4147, checked in by cetlod, 8 years ago

merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169

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