New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diadct.F90 in NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DIA/diadct.F90 @ 12097

Last change on this file since 12097 was 11536, checked in by smasson, 5 years ago

trunk: merge dev_r10984_HPC-13 into the trunk

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