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

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

changed heat calculation to use C in place of Kelvin, and initialised some previously unitialised arrays

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