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

Last change on this file since 12489 was 12489, checked in by davestorkey, 7 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

  • Property svn:keywords set to Id
File size: 58.4 KB
Line 
1MODULE diadct
2   !!======================================================================
3   !!                       ***  MODULE  diadct  ***
4   !! Ocean diagnostics: Compute the transport trough a sec.
5   !!======================================================================
6   !! History :  OPA  ! 02/1999 (Y Drillet)  original code
7   !!                 ! 10/2001 (Y Drillet, R Bourdalle Badie)
8   !!   NEMO     1.0  ! 10/2005 (M Laborie) F90
9   !!            3.0  ! 04/2007 (G Garric) Ice sections
10   !!             -   ! 04/2007 (C Bricaud) test on sec%nb_point, initialisation of ztransp1,ztransp2,...
11   !!            3.4  ! 09/2011 (C Bricaud)
12   !!----------------------------------------------------------------------
13   !! does not work with agrif
14#if ! defined key_agrif
15   !!----------------------------------------------------------------------
16   !!   dia_dct      :  Compute the transport through a sec.
17   !!   dia_dct_init :  Read namelist.
18   !!   readsec      :  Read sections description and pathway
19   !!   removepoints :  Remove points which are common to 2 procs
20   !!   transport    :  Compute transport for each sections
21   !!   dia_dct_wri  :  Write tranports results in ascii files
22   !!   interp       :  Compute temperature/salinity/density at U-point or V-point
23   !!   
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE phycst          ! physical constants
28   USE in_out_manager  ! I/O manager
29   USE daymod          ! calendar
30   USE dianam          ! build name of file
31   USE lib_mpp         ! distributed memory computing library
32#if defined key_si3
33   USE ice
34#endif
35   USE domvvl
36   USE timing          ! preformance summary
37
38   IMPLICIT NONE
39   PRIVATE
40
41   PUBLIC   dia_dct      ! routine called by step.F90
42   PUBLIC   dia_dct_init ! routine called by nemogcm.F90
43
44   !                         !!** namelist variables **
45   LOGICAL, PUBLIC ::   ln_diadct     !: Calculate transport thru a section or not
46   INTEGER         ::   nn_dct        !  Frequency of computation
47   INTEGER         ::   nn_dctwri     !  Frequency of output
48   INTEGER         ::   nn_secdebug   !  Number of the section to debug
49   
50   INTEGER, PARAMETER :: nb_class_max  = 10
51   INTEGER, PARAMETER :: nb_sec_max    = 150
52   INTEGER, PARAMETER :: nb_point_max  = 2000
53   INTEGER, PARAMETER :: nb_type_class = 10
54   INTEGER, PARAMETER :: nb_3d_vars    = 3 
55   INTEGER, PARAMETER :: nb_2d_vars    = 2 
56   INTEGER            :: nb_sec 
57
58   TYPE POINT_SECTION
59      INTEGER :: I,J
60   END TYPE POINT_SECTION
61
62   TYPE COORD_SECTION
63      REAL(wp) :: lon,lat
64   END TYPE COORD_SECTION
65
66   TYPE SECTION
67      CHARACTER(len=60)                            :: name              ! name of the sec
68      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
69                                                                       ! heat transports
70      LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
71      LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
72      TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
73      INTEGER                                      :: nb_class          ! number of boundaries for density classes
74      INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
75      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
76      REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
77                                                      zsigp           ,&! potential density classes    (99 if you don't want)
78                                                      zsal            ,&! salinity classes   (99 if you don't want)
79                                                      ztem            ,&! temperature classes(99 if you don't want)
80                                                      zlay              ! level classes      (99 if you don't want)
81      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
82      REAL(wp)                                         :: slopeSection  ! slope of the section
83      INTEGER                                          :: nb_point      ! number of points in the section
84      TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
85   END TYPE SECTION
86
87   TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
88 
89   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
90   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d 
91
92   !!----------------------------------------------------------------------
93   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
94   !! $Id$
95   !! Software governed by the CeCILL license (see ./LICENSE)
96   !!----------------------------------------------------------------------
97CONTAINS
98 
99   INTEGER FUNCTION diadct_alloc() 
100      !!----------------------------------------------------------------------
101      !!                   ***  FUNCTION diadct_alloc  ***
102      !!----------------------------------------------------------------------
103
104      ALLOCATE( transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), &
105         &      transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=diadct_alloc ) 
106
107      CALL mpp_sum( 'diadct', diadct_alloc ) 
108      IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' ) 
109
110   END FUNCTION diadct_alloc
111
112   SUBROUTINE dia_dct_init
113      !!---------------------------------------------------------------------
114      !!               ***  ROUTINE diadct  *** 
115      !!
116      !!  ** Purpose: Read the namelist parameters
117      !!              Open output files
118      !!
119      !!---------------------------------------------------------------------
120      INTEGER  ::   ios                 ! Local integer output status for namelist read
121      !!
122      NAMELIST/nam_diadct/ln_diadct, nn_dct, nn_dctwri, nn_secdebug
123      !!---------------------------------------------------------------------
124
125     READ  ( numnam_ref, nam_diadct, IOSTAT = ios, ERR = 901)
126901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadct in reference namelist' )
127
128     READ  ( numnam_cfg, nam_diadct, IOSTAT = ios, ERR = 902 )
129902  IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diadct in configuration namelist' )
130     IF(lwm) WRITE ( numond, nam_diadct )
131
132     IF( lwp ) THEN
133        WRITE(numout,*) " "
134        WRITE(numout,*) "diadct_init: compute transports through sections "
135        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
136        WRITE(numout,*) "       Calculate transport thru sections: ln_diadct = ", ln_diadct
137        WRITE(numout,*) "       Frequency of computation:          nn_dct    = ", nn_dct
138        WRITE(numout,*) "       Frequency of write:                nn_dctwri = ", nn_dctwri
139
140        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
141                                            WRITE(numout,*)"       Debug section number: ", nn_secdebug 
142        ELSE IF ( nn_secdebug ==  0 )THEN ; WRITE(numout,*)"       No section to debug"
143        ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)"       Debug all sections"
144        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
145        ENDIF
146     ENDIF
147
148     IF( ln_diadct ) THEN
149        ! control
150        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
151           &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
152
153        ! allocate dia_dct arrays
154        IF( diadct_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )
155
156        !Read section_ijglobal.diadct
157        CALL readsec
158
159        !open output file
160        IF( lwm ) THEN
161           CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
162           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
163           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
164        ENDIF
165
166        ! Initialise arrays to zero
167        transports_3d(:,:,:,:)=0.0 
168        transports_2d(:,:,:)  =0.0 
169        !
170     ENDIF
171     !
172  END SUBROUTINE dia_dct_init
173 
174 
175  SUBROUTINE dia_dct( kt, Kmm )
176     !!---------------------------------------------------------------------
177     !!               ***  ROUTINE diadct  *** 
178     !!
179     !!  Purpose :: Compute section transports and write it in numdct files
180     !!   
181     !!  Method  :: All arrays initialised to zero in dct_init
182     !!             Each nn_dct time step call subroutine 'transports' for
183     !!               each section to sum the transports over each grid cell.
184     !!             Each nn_dctwri time step:
185     !!               Divide the arrays by the number of summations to gain
186     !!               an average value
187     !!               Call dia_dct_sum to sum relevant grid boxes to obtain
188     !!               totals for each class (density, depth, temp or sal)
189     !!               Call dia_dct_wri to write the transports into file
190     !!               Reinitialise all relevant arrays to zero
191     !!---------------------------------------------------------------------
192     INTEGER, INTENT(in) ::   kt    ! ocean time step
193     INTEGER, INTENT(in) ::   Kmm   ! time level index
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(Kmm,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(Kmm,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(Kmm,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     INTEGER      ,INTENT(IN)    :: Kmm         ! time level index
580     TYPE(SECTION),INTENT(INOUT) :: sec
581     LOGICAL      ,INTENT(IN)    :: ld_debug
582     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
583     !
584     INTEGER ::   jk, jseg, jclass,jl, isgnu, isgnv    ! loop on level/segment/classes/ice categories
585     REAL(wp)::   zumid, zvmid, zumid_ice, zvmid_ice   ! U/V ocean & ice velocity on a cell segment
586     REAL(wp)::   zTnorm                               ! transport of velocity through one cell's sides
587     REAL(wp)::   ztn, zsn, zrhoi, zrhop, zsshn, zdep  ! temperature/salinity/potential density/ssh/depth at u/v point
588     TYPE(POINT_SECTION) ::   k
589      !!--------------------------------------------------------
590      !
591      IF( ld_debug )WRITE(numout,*)'      Compute transport'
592
593      !---------------------------!
594      !  COMPUTE TRANSPORT        !
595      !---------------------------!
596      IF(sec%nb_point .NE. 0)THEN   
597
598         !----------------------------------------------------------------------------------------------------
599         !Compute sign for velocities:
600         !
601         !convention:
602         !   non horizontal section: direction + is toward left hand of section
603         !       horizontal section: direction + is toward north of section
604         !
605         !
606         !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
607         !       ----------------      -----------------     ---------------             --------------
608         !
609         !   isgnv=1         direction +     
610         !  ______         _____             ______                                                   
611         !        |           //|            |                  |                         direction +   
612         !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
613         !        |_______  //         ______|    \\            | ---\                        |
614         !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
615         !               |             |          __\\|         |                   
616         !               |             |     direction +        |                      isgnv=1                                 
617         !                                                     
618         !----------------------------------------------------------------------------------------------------
619         isgnu = 1
620         IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
621         ELSE                                ; isgnv =  1
622         ENDIF
623         IF( sec%slopeSection .GE. 9999. )     isgnv =  1
624
625         IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
626
627         !--------------------------------------!
628         ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
629         !--------------------------------------!
630         DO jseg=1,MAX(sec%nb_point-1,0)
631             
632            !-------------------------------------------------------------------------------------------
633            ! Select the appropriate coordinate for computing the velocity of the segment
634            !
635            !                      CASE(0)                                    Case (2)
636            !                      -------                                    --------
637            !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
638            !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
639            !                                                                            |
640            !                                                                            |
641            !                                                                            |
642            !                      Case (3)                                            U(i,j)
643            !                      --------                                              |
644            !                                                                            |
645            !  listPoint(jseg+1) F(i,j+1)                                                |
646            !                        |                                                   |
647            !                        |                                                   |
648            !                        |                                 listPoint(jseg+1) F(i,j-1)
649            !                        |                                           
650            !                        |                                           
651            !                     U(i,j+1)                                           
652            !                        |                                       Case(1)     
653            !                        |                                       ------     
654            !                        |                                           
655            !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
656            !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
657            ! listPoint(jseg)     F(i,j)
658            !
659            !-------------------------------------------------------------------------------------------
660
661            SELECT CASE( sec%direction(jseg) )
662            CASE(0)   ;    k = sec%listPoint(jseg)
663            CASE(1)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
664            CASE(2)   ;    k = sec%listPoint(jseg)
665            CASE(3)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
666            END SELECT
667
668            !---------------------------|
669            !     LOOP ON THE LEVEL     |
670            !---------------------------|
671            DO jk = 1, mbkt(k%I,k%J)            !Sum of the transport on the vertical
672            !           ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
673            SELECT CASE( sec%direction(jseg) )
674               CASE(0,1) 
675                  ztn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) ) 
676                  zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
677                  zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
678                  zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rho0+rho0) 
679                  zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I,k%J+1,Kmm)    ) * vmask(k%I,k%J,1) 
680               CASE(2,3) 
681                  ztn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) ) 
682                  zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
683                  zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
684                  zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rho0+rho0) 
685                  zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1) 
686               END SELECT 
687               !
688               zdep= gdept(k%I,k%J,jk,Kmm) 
689 
690               SELECT CASE( sec%direction(jseg) )                !compute velocity with the correct direction
691               CASE(0,1)   
692                  zumid=0._wp
693                  zvmid=isgnv*vv(k%I,k%J,jk,Kmm)*vmask(k%I,k%J,jk) 
694               CASE(2,3) 
695                  zumid=isgnu*uu(k%I,k%J,jk,Kmm)*umask(k%I,k%J,jk) 
696                  zvmid=0._wp
697               END SELECT 
698 
699               !zTnorm=transport through one cell;
700               !velocity* cell's length * cell's thickness
701               zTnorm = zumid*e2u(k%I,k%J) * e3u(k%I,k%J,jk,Kmm)     & 
702                  &   + zvmid*e1v(k%I,k%J) * e3v(k%I,k%J,jk,Kmm) 
703
704!!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!!
705               IF( ln_linssh ) THEN              !add transport due to free surface
706                  IF( jk==1 ) THEN
707                     zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk)   & 
708                        &            + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
709                  ENDIF
710               ENDIF
711!!gm end
712              !COMPUTE TRANSPORT 
713 
714              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
715 
716              IF( sec%llstrpond ) THEN
717                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp
718                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001
719              ENDIF
720   
721           END DO !end of loop on the level
722
723#if defined key_si3
724
725           !ICE CASE   
726           !------------
727           IF( sec%ll_ice_section )THEN
728              SELECT CASE (sec%direction(jseg))
729              CASE(0)
730                 zumid_ice = 0
731                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
732              CASE(1)
733                 zumid_ice = 0
734                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
735              CASE(2)
736                 zvmid_ice = 0
737                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
738              CASE(3)
739                 zvmid_ice = 0
740                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
741              END SELECT
742   
743              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
744
745#if defined key_si3
746              DO jl=1,jpl
747                 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*       &
748                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) *  &
749                                  ( h_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  &
750                                    h_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
751                                   
752                 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
753                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
754              END DO
755#endif
756   
757           ENDIF !end of ice case
758#endif
759 
760        END DO !end of loop on the segment
761
762     ENDIF !end of sec%nb_point =0 case
763     !
764  END SUBROUTINE transport
765
766
767  SUBROUTINE dia_dct_sum(Kmm,sec,jsec) 
768     !!-------------------------------------------------------------
769     !! Purpose: Average the transport over nn_dctwri time steps 
770     !! and sum over the density/salinity/temperature/depth classes
771     !!
772     !! Method:   Sum over relevant grid cells to obtain values 
773     !!           for each class
774     !!              There are several loops:                 
775     !!              loop on the segment between 2 nodes
776     !!              loop on the level jk
777     !!              loop on the density/temperature/salinity/level classes
778     !!              test on the density/temperature/salinity/level
779     !!
780     !!  Note:    Transport through a given section is equal to the sum of transports
781     !!           computed on each proc.
782     !!           On each proc,transport is equal to the sum of transport computed through
783     !!           segments linking each point of sec%listPoint  with the next one.   
784     !!
785     !!-------------------------------------------------------------
786     INTEGER      ,INTENT(IN)    :: Kmm         ! time level index
787     TYPE(SECTION),INTENT(INOUT) :: sec 
788     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
789 
790     TYPE(POINT_SECTION) :: k 
791     INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes 
792     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point
793     !!-------------------------------------------------------------
794 
795     !! Sum the relevant segments to obtain values for each class
796     IF(sec%nb_point .NE. 0)THEN   
797 
798        !--------------------------------------!
799        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
800        !--------------------------------------!
801        DO jseg=1,MAX(sec%nb_point-1,0) 
802           
803           !-------------------------------------------------------------------------------------------
804           ! Select the appropriate coordinate for computing the velocity of the segment
805           !
806           !                      CASE(0)                                    Case (2)
807           !                      -------                                    --------
808           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
809           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
810           !                                                                            |
811           !                                                                            |
812           !                                                                            |
813           !                      Case (3)                                            U(i,j)
814           !                      --------                                              |
815           !                                                                            |
816           !  listPoint(jseg+1) F(i,j+1)                                                |
817           !                        |                                                   |
818           !                        |                                                   |
819           !                        |                                 listPoint(jseg+1) F(i,j-1)
820           !                        |                                             
821           !                        |                                             
822           !                     U(i,j+1)                                             
823           !                        |                                       Case(1)     
824           !                        |                                       ------       
825           !                        |                                             
826           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
827           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
828           ! listPoint(jseg)     F(i,j)
829           
830           !-------------------------------------------------------------------------------------------
831 
832           SELECT CASE( sec%direction(jseg) ) 
833           CASE(0)  ;   k = sec%listPoint(jseg) 
834           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
835           CASE(2)  ;   k = sec%listPoint(jseg) 
836           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
837           END SELECT 
838 
839           !---------------------------|
840           !     LOOP ON THE LEVEL     |
841           !---------------------------|
842           !Sum of the transport on the vertical 
843           DO jk=1,mbkt(k%I,k%J) 
844 
845              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
846              SELECT CASE( sec%direction(jseg) ) 
847              CASE(0,1) 
848                 ztn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) ) 
849                 zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
850                 zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
851                 zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rho0+rho0) 
852
853              CASE(2,3) 
854                 ztn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) ) 
855                 zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
856                 zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
857                 zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rho0+rho0) 
858                 zsshn =  0.5*( ssh(k%I,k%J,Kmm)    + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1) 
859              END SELECT
860 
861              zdep= gdept(k%I,k%J,jk,Kmm) 
862 
863              !-------------------------------
864              !  LOOP ON THE DENSITY CLASSES |
865              !-------------------------------
866              !The computation is made for each density/temperature/salinity/depth class
867              DO jclass=1,MAX(1,sec%nb_class-1) 
868 
869                 !----------------------------------------------!
870                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 
871                 !----------------------------------------------!
872
873                 IF ( (                                                    & 
874                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
875                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
876                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
877 
878                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
879                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
880                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
881 
882                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
883                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
884                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
885 
886                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
887                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
888                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
889 
890                    ((( zdep .GE. sec%zlay(jclass)) .AND.                & 
891                    (   zdep .LE. sec%zlay(jclass+1))) .OR.              & 
892                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
893                                                                   ))   THEN 
894 
895                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
896                    !----------------------------------------------------------------------------
897                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN 
898                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
899                    ELSE
900                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
901                    ENDIF
902                    IF( sec%llstrpond )THEN
903 
904                       IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN
905                          sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk) 
906                       ELSE
907                          sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk) 
908                       ENDIF
909 
910                       IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN
911                          sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk) 
912                       ELSE
913                          sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk) 
914                       ENDIF
915 
916                    ELSE
917                       sec%transport( 3,jclass) = 0._wp 
918                       sec%transport( 4,jclass) = 0._wp 
919                       sec%transport( 5,jclass) = 0._wp 
920                       sec%transport( 6,jclass) = 0._wp 
921                    ENDIF
922 
923                 ENDIF ! end of test if point is in class
924   
925              END DO ! end of loop on the classes
926 
927           END DO ! loop over jk
928 
929#if defined key_si3 
930 
931           !ICE CASE     
932           IF( sec%ll_ice_section )THEN
933 
934              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
935                 sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6 
936              ELSE
937                 sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6 
938              ENDIF
939 
940              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
941                 sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6 
942              ELSE
943                 sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6 
944              ENDIF
945 
946           ENDIF !end of ice case
947#endif
948 
949        END DO !end of loop on the segment
950 
951     ELSE  !if sec%nb_point =0
952        sec%transport(1:2,:)=0. 
953        IF (sec%llstrpond) sec%transport(3:6,:)=0. 
954        IF (sec%ll_ice_section) sec%transport(7:10,:)=0. 
955     ENDIF !end of sec%nb_point =0 case
956 
957  END SUBROUTINE dia_dct_sum 
958
959
960  SUBROUTINE dia_dct_wri(kt,ksec,sec)
961     !!-------------------------------------------------------------
962     !! Write transport output in numdct
963     !!
964     !! Purpose: Write  transports in ascii files
965     !!
966     !! Method:
967     !!        1. Write volume transports in "volume_transport"
968     !!           Unit: Sv : area * Velocity / 1.e6
969     !!
970     !!        2. Write heat transports in "heat_transport"
971     !!           Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
972     !!
973     !!        3. Write salt transports in "salt_transport"
974     !!           Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9
975     !!
976     !!-------------------------------------------------------------
977     !!arguments
978     INTEGER, INTENT(IN)          :: kt          ! time-step
979     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
980     INTEGER ,INTENT(IN)          :: ksec        ! section number
981
982     !!local declarations
983     INTEGER               :: jclass             ! Dummy loop
984     CHARACTER(len=2)      :: classe             ! Classname
985     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
986     REAL(wp)              :: zslope             ! section's slope coeff
987     !
988     REAL(wp), DIMENSION(nb_type_class)::   zsumclasses   ! 1D workspace
989     !!-------------------------------------------------------------
990
991     zsumclasses(:)=0._wp
992     zslope = sec%slopeSection       
993
994 
995     DO jclass=1,MAX(1,sec%nb_class-1)
996
997        classe   = 'N       '
998        zbnd1   = 0._wp
999        zbnd2   = 0._wp
1000        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
1001
1002   
1003        !insitu density classes transports
1004        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
1005            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
1006           classe = 'DI       '
1007           zbnd1 = sec%zsigi(jclass)
1008           zbnd2 = sec%zsigi(jclass+1)
1009        ENDIF
1010        !potential density classes transports
1011        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
1012            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
1013           classe = 'DP      '
1014           zbnd1 = sec%zsigp(jclass)
1015           zbnd2 = sec%zsigp(jclass+1)
1016        ENDIF
1017        !depth classes transports
1018        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
1019            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
1020           classe = 'Z       '
1021           zbnd1 = sec%zlay(jclass)
1022           zbnd2 = sec%zlay(jclass+1)
1023        ENDIF
1024        !salinity classes transports
1025        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
1026            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
1027           classe = 'S       '
1028           zbnd1 = sec%zsal(jclass)
1029           zbnd2 = sec%zsal(jclass+1)   
1030        ENDIF
1031        !temperature classes transports
1032        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
1033            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
1034           classe = 'T       '
1035           zbnd1 = sec%ztem(jclass)
1036           zbnd2 = sec%ztem(jclass+1)
1037        ENDIF
1038                 
1039        !write volume transport per class
1040        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1041                              jclass,classe,zbnd1,zbnd2,&
1042                              sec%transport(1,jclass),sec%transport(2,jclass), &
1043                              sec%transport(1,jclass)+sec%transport(2,jclass)
1044
1045        IF( sec%llstrpond )THEN
1046
1047           !write heat transport per class:
1048           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  &
1049                              jclass,classe,zbnd1,zbnd2,&
1050                              sec%transport(3,jclass)*1.e-15,sec%transport(4,jclass)*1.e-15, &
1051                              ( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15
1052           !write salt transport per class
1053           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  &
1054                              jclass,classe,zbnd1,zbnd2,&
1055                              sec%transport(5,jclass)*1.e-9,sec%transport(6,jclass)*1.e-9,&
1056                              (sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9
1057        ENDIF
1058
1059     ENDDO
1060
1061     zbnd1 = 0._wp
1062     zbnd2 = 0._wp
1063     jclass=0
1064
1065     !write total volume transport
1066     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1067                           jclass,"total",zbnd1,zbnd2,&
1068                           zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
1069
1070     IF( sec%llstrpond )THEN
1071
1072        !write total heat transport
1073        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1074                           jclass,"total",zbnd1,zbnd2,&
1075                           zsumclasses(3)*1.e-15,zsumclasses(4)*1.e-15,&
1076                           (zsumclasses(3)+zsumclasses(4) )*1.e-15
1077        !write total salt transport
1078        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1079                           jclass,"total",zbnd1,zbnd2,&
1080                           zsumclasses(5)*1.e-9,zsumclasses(6)*1.e-9,&
1081                           (zsumclasses(5)+zsumclasses(6))*1.e-9
1082     ENDIF
1083
1084     
1085     IF ( sec%ll_ice_section) THEN
1086        !write total ice volume transport
1087        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1088                              jclass,"ice_vol",zbnd1,zbnd2,&
1089                              sec%transport(7,1),sec%transport(8,1),&
1090                              sec%transport(7,1)+sec%transport(8,1)
1091        !write total ice surface transport
1092        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1093                              jclass,"ice_surf",zbnd1,zbnd2,&
1094                              sec%transport(9,1),sec%transport(10,1), &
1095                              sec%transport(9,1)+sec%transport(10,1) 
1096      ENDIF
1097                                             
1098118   FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
1099119   FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
1100      !
1101   END SUBROUTINE dia_dct_wri
1102
1103
1104   FUNCTION interp(Kmm, ki, kj, kk, cd_point, ptab)
1105  !!----------------------------------------------------------------------
1106  !!
1107  !!   Purpose: compute temperature/salinity/density at U-point or V-point
1108  !!   --------
1109  !!
1110  !!   Method:
1111  !!   ------
1112  !!
1113  !!   ====> full step and partial step
1114  !!
1115  !!
1116  !!    |    I          |    I+1           |    Z=temperature/salinity/density at U-poinT
1117  !!    |               |                  |
1118  !!  ----------------------------------------  1. Veritcal interpolation: compute zbis
1119  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
1120  !!    |               |                  |       zbis =
1121  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
1122  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
1123  !!    |               |                  |
1124  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
1125  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K) 
1126  !!    |     .         |                  |
1127  !!    |     .         |                  |       interp = ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)
1128  !!    |     .         |                  |
1129  !!  ------------------------------------------
1130  !!    |     .         |                  |
1131  !!    |     .         |                  |
1132  !!    |     .         |                  |
1133  !!K   |    zbis.......U...ptab(I+1,J,K)  |
1134  !!    |     .         |                  |
1135  !!    | ptab(I,J,K)   |                  |
1136  !!    |               |------------------|
1137  !!    |               | partials         |
1138  !!    |               |  steps           |
1139  !!  -------------------------------------------
1140  !!    <----zet1------><----zet2--------->
1141  !!
1142  !!
1143  !!   ====>  s-coordinate
1144  !!     
1145  !!    |                |                  |   1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
1146  !!    |                |                  |      Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
1147  !!    |                | ptab(I+1,J,K)    |
1148  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point
1149  !!    |                |      ^           |   
1150  !!    |                |      | zdep2     |   
1151  !!    |                |      |           |   
1152  !!    |       ^        U      v           |
1153  !!    |       |        |                  |
1154  !!    |       | zdep1  |                  |   
1155  !!    |       v        |                  |
1156  !!    |      T1        |                  |
1157  !!    | ptab(I,J,K)    |                  |
1158  !!    |                |                  |
1159  !!    |                |                  |
1160  !!
1161  !!    <----zet1--------><----zet2--------->
1162  !!
1163  !!----------------------------------------------------------------------
1164  !*arguments
1165  INTEGER, INTENT(IN)                          :: Kmm          ! time level index
1166  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
1167  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
1168  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk )
1169  REAL(wp)                                     :: interp       ! interpolated variable
1170
1171  !*local declations
1172  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
1173  REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu            ! local real
1174  REAL(wp):: zet1, zet2                                        ! weight for interpolation
1175  REAL(wp):: zdep1,zdep2                                       ! differences of depth
1176  REAL(wp):: zmsk                                              ! mask value
1177  !!----------------------------------------------------------------------
1178
1179  IF( cd_point=='U' )THEN
1180     ii1 = ki    ; ij1 = kj 
1181     ii2 = ki+1  ; ij2 = kj 
1182
1183     zet1=e1t(ii1,ij1)
1184     zet2=e1t(ii2,ij2)
1185     zmsk=umask(ii1,ij1,kk)
1186 
1187
1188  ELSE ! cd_point=='V'
1189     ii1 = ki    ; ij1 = kj 
1190     ii2 = ki    ; ij2 = kj+1 
1191
1192     zet1=e2t(ii1,ij1)
1193     zet2=e2t(ii2,ij2)
1194     zmsk=vmask(ii1,ij1,kk)
1195
1196  ENDIF
1197
1198  IF( ln_sco )THEN   ! s-coordinate case
1199
1200     zdepu = ( gdept(ii1,ij1,kk,Kmm) +  gdept(ii2,ij2,kk,Kmm) ) * 0.5_wp 
1201     zdep1 = gdept(ii1,ij1,kk,Kmm) - zdepu
1202     zdep2 = gdept(ii2,ij2,kk,Kmm) - zdepu
1203
1204     ! weights
1205     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1206     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1207 
1208     ! result
1209     interp = zmsk * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )   
1210
1211
1212  ELSE       ! full step or partial step case
1213
1214     ze3t  = e3t(ii2,ij2,kk,Kmm) - e3t(ii1,ij1,kk,Kmm) 
1215     zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) ) / e3w(ii2,ij2,kk,Kmm)
1216     zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) ) / e3w(ii1,ij1,kk,Kmm)
1217
1218     IF(kk .NE. 1)THEN
1219
1220        IF( ze3t >= 0. )THEN 
1221           ! zbis
1222           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
1223           ! result
1224            interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
1225        ELSE
1226           ! zbis
1227           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1228           ! result
1229           interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1230        ENDIF   
1231
1232     ELSE
1233        interp = zmsk * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1234     ENDIF
1235
1236  ENDIF
1237      !
1238   END FUNCTION interp
1239
1240#else
1241   !!----------------------------------------------------------------------
1242   !!   Dummy module                                             
1243   !!----------------------------------------------------------------------
1244   LOGICAL, PUBLIC ::   ln_diadct = .FALSE.
1245CONTAINS
1246   SUBROUTINE dia_dct_init
1247      IMPLICIT NONE
1248   END SUBROUTINE dia_dct_init
1249
1250   SUBROUTINE dia_dct( kt, Kmm )         ! Dummy routine
1251      IMPLICIT NONE
1252      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
1253      INTEGER, INTENT( in ) :: Kmm  ! ocean time level index
1254      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
1255   END SUBROUTINE dia_dct
1256   !
1257#endif
1258
1259   !!======================================================================
1260END MODULE diadct
Note: See TracBrowser for help on using the repository browser.