source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_DYN_optimization/src/OCE/DIA/diadct.F90 @ 11381

Last change on this file since 11381 was 11381, checked in by girrmann, 2 years ago

dev_r10984_HPC-13 : repair faulty commit [11380], see #2308

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