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

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

diadct: correct bugs on directions. ticket 986

  • Property svn:executable set to *
File size: 53.7 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 .GE. 0 ) THEN  ; isgnv = -1 
625        ELSE                                ; isgnv =  1
626        ENDIF
627
628        IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
629
630        !--------------------------------------!
631        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
632        !--------------------------------------!
633        DO jseg=1,MAX(sec%nb_point-1,0)
634             
635           !-------------------------------------------------------------------------------------------
636           ! Select the appropriate coordinate for computing the velocity of the segment
637           !
638           !                      CASE(0)                                    Case (2)
639           !                      -------                                    --------
640           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
641           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
642           !                                                                            |
643           !                                                                            |
644           !                                                                            |
645           !                      Case (3)                                            U(i,j)
646           !                      --------                                              |
647           !                                                                            |
648           !  listPoint(jseg+1) F(i,j+1)                                                |
649           !                        |                                                   |
650           !                        |                                                   |
651           !                        |                                 listPoint(jseg+1) F(i,j-1)
652           !                        |                                           
653           !                        |                                           
654           !                     U(i,j+1)                                           
655           !                        |                                       Case(1)     
656           !                        |                                       ------     
657           !                        |                                           
658           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
659           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
660           ! listPoint(jseg)     F(i,j)
661           !
662           !-------------------------------------------------------------------------------------------
663
664           SELECT CASE( sec%direction(jseg) )
665           CASE(0)  ;   k = sec%listPoint(jseg)
666           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
667           CASE(2)  ;   k = sec%listPoint(jseg)
668           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
669           END SELECT
670
671           !-------------------------------
672           !  LOOP ON THE DENSITY CLASSES |
673           !-------------------------------
674           !The computation is made for each density class
675           DO jclass=1,MAX(1,sec%nb_class-1)
676
677              ztransp1=0._wp ; zrhoi1=0._wp ; zrhop1=0._wp ; ztemp1=0._wp ;zsal1=0._wp
678              ztransp2=0._wp ; zrhoi2=0._wp ; zrhop2=0._wp ; ztemp2=0._wp ;zsal2=0._wp
679   
680              !---------------------------|
681              !     LOOP ON THE LEVEL     |
682              !---------------------------|
683              !Sum of the transport on the vertical
684              DO jk=1,jpk
685                   
686
687                 ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point
688                 SELECT CASE( sec%direction(jseg) )
689                 CASE(0,1)
690                    ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
691                    zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
692                    zrhop = interp(k%I,k%J,jk,'V',rhop)
693                    zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
694                    zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
695                 CASE(2,3)
696                    ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
697                    zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
698                    zrhop = interp(k%I,k%J,jk,'U',rhop)
699                    zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
700                    zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
701                 END SELECT
702
703                 zfsdep= gdept(k%I,k%J,jk)
704 
705                 !----------------------------------------------!
706                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
707                 !----------------------------------------------!
708 
709                 IF ( (    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.    &
710                           (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.    &
711                           ( sec%zsigp(jclass) .EQ. 99.)) .AND.                 &
712                           ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
713                           (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.    &
714                           ( sec%zsigi(jclass) .EQ. 99.)) .AND.                 &
715                           ((( zsn .GT. sec%zsal(jclass)) .AND.                &
716                           (   zsn .LE. sec%zsal(jclass+1))) .OR.              &
717                           ( sec%zsal(jclass) .EQ. 99.)) .AND.                 &
718                           ((( ztn .GE. sec%ztem(jclass)) .AND.                &
719                           (   ztn .LE. sec%ztem(jclass+1))) .OR.              &
720                           ( sec%ztem(jclass) .EQ.99.)) .AND.                  &
721                           ((( zfsdep .GE. sec%zlay(jclass)) .AND.            &
722                           (   zfsdep .LE. sec%zlay(jclass+1))) .OR.          &
723                           ( sec%zlay(jclass) .EQ. 99. ))))   THEN
724
725
726                    !compute velocity with the correct direction
727                    SELECT CASE( sec%direction(jseg) )
728                    CASE(0,1) 
729                       zumid=0.
730                       zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
731                    CASE(2,3)
732                       zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
733                       zvmid=0.
734                    END SELECT
735
736                    !velocity* cell's length * cell's thickness
737                    zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &
738                           zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)
739
740#if ! defined key_vvl
741                    !add transport due to free surface
742                    IF( jk==1 )THEN
743                       zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
744                                         zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
745                    ENDIF
746#endif
747                    !COMPUTE TRANSPORT
748                    !zTnorm=transport through one cell for one class
749                    !ztransp1 or ztransp2=transport through one cell i
750                    !                     for one class for one direction
751                    IF( zTnorm .GE. 0 )THEN
752
753                       ztransp1=zTnorm+ztransp1
754 
755                       IF ( sec%llstrpond ) THEN
756                          ztemp1 = ztemp1  + zTnorm * ztn 
757                          zsal1  = zsal1   + zTnorm * zsn
758                          zrhoi1 = zrhoi1  + zTnorm * zrhoi
759                          zrhop1 = zrhop1  + zTnorm * zrhop
760                       ENDIF
761
762                    ELSE
763
764                       ztransp2=(zTnorm)+ztransp2
765
766                       IF ( sec%llstrpond ) THEN
767                          ztemp2 = ztemp2  + zTnorm * ztn 
768                          zsal2  = zsal2   + zTnorm * zsn
769                          zrhoi2 = zrhoi2  + zTnorm * zrhoi
770                          zrhop2 = zrhop2  + zTnorm * zrhop
771                       ENDIF
772                    ENDIF
773 
774           
775                 ENDIF ! end of density test
776              ENDDO!end of loop on the level
777
778              !ZSUM=TRANSPORT FOR EACH CLASSES FOR THE  DIRECTIONS
779              !---------------------------------------------------
780              zsum(1,jclass)     = zsum(1,jclass)+ztransp1
781              zsum(2,jclass)     = zsum(2,jclass)+ztransp2
782              IF( sec%llstrpond )THEN
783                 zsum(3 ,jclass) = zsum( 3,jclass)+zrhoi1
784                 zsum(4 ,jclass) = zsum( 4,jclass)+zrhoi2
785                 zsum(5 ,jclass) = zsum( 5,jclass)+zrhop1
786                 zsum(6 ,jclass) = zsum( 6,jclass)+zrhop2
787                 zsum(7 ,jclass) = zsum( 7,jclass)+ztemp1
788                 zsum(8 ,jclass) = zsum( 8,jclass)+ztemp2
789                 zsum(9 ,jclass) = zsum( 9,jclass)+zsal1
790                 zsum(10,jclass) = zsum(10,jclass)+zsal2
791              ENDIF
792   
793           ENDDO !end of loop on the density classes
794
795#if defined key_lim2 || defined key_lim3
796
797           !ICE CASE   
798           !------------
799           IF( sec%ll_ice_section )THEN
800              SELECT CASE (sec%direction(jseg))
801              CASE(0)
802                 zumid_ice = 0
803                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
804              CASE(1)
805                 zumid_ice = 0
806                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
807              CASE(2)
808                 zvmid_ice = 0
809                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
810              CASE(3)
811                 zvmid_ice = 0
812                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
813              END SELECT
814   
815              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
816   
817              IF( zTnorm .GE. 0)THEN
818                 zice_vol_pos = (zTnorm)*   &
819                                      (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
820                                     *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
821                                       hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
822                                      +zice_vol_pos
823                 zice_surf_pos = (zTnorm)*   &
824                                       (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
825                                      +zice_surf_pos
826              ELSE
827                 zice_vol_neg=(zTnorm)*   &
828                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
829                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
830                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
831                                  +zice_vol_neg
832                 zice_surf_neg=(zTnorm)*   &
833                                    (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
834                                     +zice_surf_neg
835              ENDIF
836   
837              zsum(11,1) = zsum(11,1)+zice_vol_pos
838              zsum(12,1) = zsum(12,1)+zice_vol_neg
839              zsum(13,1) = zsum(13,1)+zice_surf_pos
840              zsum(14,1) = zsum(14,1)+zice_surf_neg
841   
842           ENDIF !end of ice case
843#endif
844 
845        ENDDO !end of loop on the segment
846
847
848     ELSE  !if sec%nb_point =0
849        zsum(1:2,:)=0.
850        IF (sec%llstrpond) zsum(3:10,:)=0.
851        zsum( 11:14,:)=0.
852     ENDIF   !end of sec%nb_point =0 case
853
854     !-------------------------------|
855     !FINISH COMPUTING TRANSPORTS    |
856     !-------------------------------|
857     DO jclass=1,MAX(1,sec%nb_class-1)
858        sec%transport(1,jclass)=sec%transport(1,jclass)+zsum(1,jclass)*1.E-6
859        sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6
860        IF( sec%llstrpond ) THEN
861           IF( zsum(1,jclass) .NE. 0._wp ) THEN
862              sec%transport( 3,jclass) = sec%transport( 3,jclass) + zsum( 3,jclass)/zsum(1,jclass)
863              sec%transport( 5,jclass) = sec%transport( 5,jclass) + zsum( 5,jclass)/zsum(1,jclass)
864              sec%transport( 7,jclass) = sec%transport( 7,jclass) + zsum( 7,jclass)
865              sec%transport( 9,jclass) = sec%transport( 9,jclass) + zsum( 9,jclass)
866           ENDIF
867           IF( zsum(2,jclass) .NE. 0._wp )THEN
868              sec%transport( 4,jclass) = sec%transport( 4,jclass) + zsum( 4,jclass)/zsum(2,jclass)
869              sec%transport( 6,jclass) = sec%transport( 6,jclass) + zsum( 6,jclass)/zsum(2,jclass)
870              sec%transport( 8,jclass) = sec%transport( 8,jclass) + zsum( 8,jclass)
871              sec%transport(10,jclass) = sec%transport(10,jclass) + zsum(10,jclass)
872           ENDIF
873        ELSE
874           sec%transport( 3,jclass) = 0._wp
875           sec%transport( 4,jclass) = 0._wp
876           sec%transport( 5,jclass) = 0._wp
877           sec%transport( 6,jclass) = 0._wp
878           sec%transport( 7,jclass) = 0._wp
879           sec%transport( 8,jclass) = 0._wp
880           sec%transport(10,jclass) = 0._wp
881        ENDIF
882     ENDDO   
883
884     IF( sec%ll_ice_section ) THEN
885        sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6
886        sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6
887        sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6
888        sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6
889     ENDIF
890
891     CALL wrk_dealloc( nb_type_class , nb_class_max , zsum   )
892     !
893  END SUBROUTINE transport
894 
895  SUBROUTINE dia_dct_wri(kt,ksec,sec)
896     !!-------------------------------------------------------------
897     !! Write transport output in numdct
898     !!
899     !! Purpose: Write  transports in ascii files
900     !!
901     !! Method:
902     !!        1. Write volume transports in "volume_transport"
903     !!           Unit: Sv : area * Velocity / 1.e6
904     !!
905     !!        2. Write heat transports in "heat_transport"
906     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
907     !!
908     !!        3. Write salt transports in "salt_transport"
909     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
910     !!
911     !!-------------------------------------------------------------
912     !!arguments
913     INTEGER, INTENT(IN)          :: kt          ! time-step
914     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
915     INTEGER ,INTENT(IN)          :: ksec        ! section number
916
917     !!local declarations
918     INTEGER               :: jcl,ji             ! Dummy loop
919     CHARACTER(len=2)      :: classe             ! Classname
920     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
921     REAL(wp)              :: zslope             ! section's slope coeff
922     !
923     REAL(wp), POINTER, DIMENSION(:):: zsumclass ! 1D workspace
924     !!-------------------------------------------------------------
925     CALL wrk_alloc(nb_type_class , zsumclass ) 
926
927     zsumclass(:)=0._wp
928     zslope = sec%slopeSection       
929
930 
931     DO jcl=1,MAX(1,sec%nb_class-1)
932
933        ! Mean computation
934        sec%transport(:,jcl)=sec%transport(:,jcl)/(nn_dctwri/nn_dct)
935        classe   = 'N       '
936        zbnd1   = 0._wp
937        zbnd2   = 0._wp
938        zsumclass(1:nb_type_class)=zsumclass(1:nb_type_class)+sec%transport(1:nb_type_class,jcl)
939
940   
941        !insitu density classes transports
942        IF( ( sec%zsigi(jcl)   .NE. 99._wp ) .AND. &
943            ( sec%zsigi(jcl+1) .NE. 99._wp )       )THEN
944           classe = 'DI       '
945           zbnd1 = sec%zsigi(jcl)
946           zbnd2 = sec%zsigi(jcl+1)
947        ENDIF
948        !potential density classes transports
949        IF( ( sec%zsigp(jcl)   .NE. 99._wp ) .AND. &
950            ( sec%zsigp(jcl+1) .NE. 99._wp )       )THEN
951           classe = 'DP      '
952           zbnd1 = sec%zsigp(jcl)
953           zbnd2 = sec%zsigp(jcl+1)
954        ENDIF
955        !depth classes transports
956        IF( ( sec%zlay(jcl)    .NE. 99._wp ) .AND. &
957            ( sec%zlay(jcl+1)  .NE. 99._wp )       )THEN
958           classe = 'Z       '
959           zbnd1 = sec%zlay(jcl)
960           zbnd2 = sec%zlay(jcl+1)
961        ENDIF
962        !salinity classes transports
963        IF( ( sec%zsal(jcl) .NE. 99._wp    ) .AND. &
964            ( sec%zsal(jcl+1) .NE. 99._wp  )       )THEN
965           classe = 'S       '
966           zbnd1 = sec%zsal(jcl)
967           zbnd2 = sec%zsal(jcl+1)   
968        ENDIF
969        !temperature classes transports
970        IF( ( sec%ztem(jcl) .NE. 99._wp     ) .AND. &
971            ( sec%ztem(jcl+1) .NE. 99._wp     )       ) THEN
972           classe = 'T       '
973           zbnd1 = sec%ztem(jcl)
974           zbnd2 = sec%ztem(jcl+1)
975        ENDIF
976                 
977        !write volume transport per class
978        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
979                              jcl,classe,zbnd1,zbnd2,&
980                              sec%transport(1,jcl),sec%transport(2,jcl), &
981                              sec%transport(1,jcl)+sec%transport(2,jcl)
982
983        IF( sec%llstrpond )THEN
984
985           !write heat transport per class:
986           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  &
987                              jcl,classe,zbnd1,zbnd2,&
988                              sec%transport(7,jcl)*1000._wp*rcp/1.e15,sec%transport(8,jcl)*1000._wp*rcp/1.e15, &
989                              ( sec%transport(7,jcl)+sec%transport(8,jcl) )*1000._wp*rcp/1.e15
990           !write salt transport per class
991           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  &
992                              jcl,classe,zbnd1,zbnd2,&
993                              sec%transport(9,jcl)*1000._wp/1.e9,sec%transport(10,jcl)*1000._wp/1.e9,&
994                              (sec%transport(9,jcl)+sec%transport(10,jcl))*1000._wp/1.e9
995        ENDIF
996
997     ENDDO
998
999     zbnd1 = 0._wp
1000     zbnd2 = 0._wp
1001     jcl=0
1002
1003     !write total volume transport
1004     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1005                           jcl,"total",zbnd1,zbnd2,&
1006                           zsumclass(1),zsumclass(2),zsumclass(1)+zsumclass(2)
1007
1008     IF( sec%llstrpond )THEN
1009
1010        !write total heat transport
1011        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1012                           jcl,"total",zbnd1,zbnd2,&
1013                           zsumclass(7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,&
1014                           (zsumclass(7)+zsumclass(8) )* 1000._wp*rcp/1.e15
1015        !write total salt transport
1016        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1017                           jcl,"total",zbnd1,zbnd2,&
1018                           zsumclass(9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,&
1019                           (zsumclass(9)+zsumclass(10))*1000._wp/1.e9
1020     ENDIF
1021
1022     
1023     IF ( sec%ll_ice_section) THEN
1024        !write total ice volume transport
1025        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1026                              jcl,"ice_vol",zbnd1,zbnd2,&
1027                              sec%transport(9,1),sec%transport(10,1),&
1028                              sec%transport(9,1)+sec%transport(10,1)
1029        !write total ice surface transport
1030        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1031                              jcl,"ice_surf",zbnd1,zbnd2,&
1032                              sec%transport(11,1),sec%transport(12,1), &
1033                              sec%transport(11,1)+sec%transport(12,1) 
1034     ENDIF
1035                                             
1036118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
1037119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
1038
1039     CALL wrk_dealloc(nb_type_class , zsumclass ) 
1040  END SUBROUTINE dia_dct_wri
1041
1042  FUNCTION interp(ki, kj, kk, cd_point, ptab)
1043  !!----------------------------------------------------------------------
1044  !!
1045  !!   Purpose: compute Temperature/Salinity/density at U-point or V-point
1046  !!   --------
1047  !!
1048  !!   Method:
1049  !!   ------
1050  !!
1051  !!   ====> full step and partial step
1052  !!
1053  !!
1054  !!    |    I          |    I+1           |    Z=Temperature/Salinity/density at U-poinT
1055  !!    |               |                  |
1056  !!  ----------------------------------------  1. Veritcale interpolation: compute zbis
1057  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
1058  !!    |               |                  |       zbis =
1059  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
1060  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
1061  !!    |               |                  |
1062  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
1063  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K) 
1064  !!    |     .         |                  |
1065  !!    |     .         |                  |       interp = ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)
1066  !!    |     .         |                  |
1067  !!  ------------------------------------------
1068  !!    |     .         |                  |
1069  !!    |     .         |                  |
1070  !!    |     .         |                  |
1071  !!K   |    zbis.......U...ptab(I+1,J,K)  |
1072  !!    |     .         |                  |
1073  !!    | ptab(I,J,K)   |                  |
1074  !!    |               |------------------|
1075  !!    |               | partials         |
1076  !!    |               |  steps           |
1077  !!  -------------------------------------------
1078  !!    <----zet1------><----zet2--------->
1079  !!
1080  !!
1081  !!   ====>  s-coordinate
1082  !!     
1083  !!    |                |                  |   1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
1084  !!    |                |                  |      Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
1085  !!    |                | ptab(I+1,J,K)    |
1086  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point
1087  !!    |                |      ^           |   
1088  !!    |                |      | zdep2     |   
1089  !!    |                |      |           |   
1090  !!    |       ^        U      v           |
1091  !!    |       |        |                  |
1092  !!    |       | zdep1  |                  |   
1093  !!    |       v        |                  |
1094  !!    |      T1        |                  |
1095  !!    | ptab(I,J,K)    |                  |
1096  !!    |                |                  |
1097  !!    |                |                  |
1098  !!
1099  !!    <----zet1--------><----zet2--------->
1100  !!
1101  !!----------------------------------------------------------------------
1102  !*arguments
1103  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
1104  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
1105  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk )
1106  REAL(wp)                                     :: interp       ! interpolated variable
1107
1108  !*local declations
1109  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
1110  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real
1111  REAL(wp):: zet1, zet2                                        ! weight for interpolation
1112  REAL(wp):: zdep1,zdep2                                       ! differences of depth
1113  !!----------------------------------------------------------------------
1114
1115  IF( cd_point=='U' )THEN
1116     ii1 = ki    ; ij1 = kj 
1117     ii2 = ki+1  ; ij2 = kj 
1118
1119     zet1=e1t(ii1,ij1)
1120     zet2=e1t(ii2,ij2)
1121
1122  ELSE ! cd_point=='V'
1123     ii1 = ki    ; ij1 = kj 
1124     ii2 = ki    ; ij2 = kj+1 
1125
1126     zet1=e2t(ii1,ij1)
1127     zet2=e2t(ii2,ij2)
1128
1129  ENDIF
1130
1131  IF( ln_sco )THEN   ! s-coordinate case
1132
1133     zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2 
1134     zdep1 = fsdept(ii1,ij1,kk) - zdepu
1135     zdep2 = fsdept(ii2,ij2,kk) - zdepu
1136
1137     !weights
1138     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1139     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1140 
1141     ! result
1142     interp = umask(ii1,ij1,kk) * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )   
1143
1144
1145  ELSE       ! full step or partial step case
1146
1147#if defined key_vvl
1148
1149     ze3t  = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk) 
1150     zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk)
1151     zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk)
1152
1153#else
1154
1155     ze3t  = fse3t(ii2,ij2,kk)   - fse3t(ii1,ij1,kk) 
1156     zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
1157     zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
1158
1159#endif
1160
1161     IF(kk .NE. 1)THEN
1162
1163        IF( ze3t >= 0. )THEN 
1164           !zbis
1165           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
1166           ! result
1167            interp = umask(ii1,ij1,kk) * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
1168        ELSE
1169           !zbis
1170           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1171           ! result
1172           interp = umask(ii1,ij1,kk) * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1173        ENDIF   
1174
1175     ELSE
1176        interp = umask(ii1,ij1,kk) * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1177     ENDIF
1178
1179  ENDIF
1180
1181
1182  END FUNCTION interp
1183
1184#else
1185   !!----------------------------------------------------------------------
1186   !!   Default option :                                       Dummy module
1187   !!----------------------------------------------------------------------
1188   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
1189   PUBLIC
1190CONTAINS
1191
1192   SUBROUTINE dia_dct_init          ! Dummy routine
1193      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
1194   END SUBROUTINE dia_dct_init
1195
1196   SUBROUTINE dia_dct( kt )           ! Dummy routine
1197      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
1198      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
1199   END SUBROUTINE dia_dct
1200#endif
1201
1202END MODULE diadct
Note: See TracBrowser for help on using the repository browser.