source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

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