source: branches/2012/dev_r3452_UKMO2_DIADCT/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 3522

Last change on this file since 3522 was 3522, checked in by rfurner, 8 years ago

removed unneccessary calculation of mass transports and corrected small ice bug

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