source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 3452

Last change on this file since 3452 was 3452, checked in by cbricaud, 9 years ago

correction of sign for horizontal section

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