source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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