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/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DIA – NEMO

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/DIA/diadct.F90 @ 10297

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

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of mppmin/max/sum, see #2133

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