source: NEMO/trunk/src/OCE/DIA/diadct.F90 @ 13237

Last change on this file since 13237 was 13237, checked in by smasson, 3 months ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

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