source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diadct.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh

  • Property svn:keywords set to Id
File size: 58.1 KB
Line 
1MODULE diadct
2   !!======================================================================
3   !!                       ***  MODULE  diadct  ***
4   !! Ocean diagnostics: Compute the transport trough a sec.
5   !!======================================================================
6   !! History :  OPA  ! 02/1999 (Y Drillet)  original code
7   !!                 ! 10/2001 (Y Drillet, R Bourdalle Badie)
8   !!   NEMO     1.0  ! 10/2005 (M Laborie) F90
9   !!            3.0  ! 04/2007 (G Garric) Ice sections
10   !!             -   ! 04/2007 (C Bricaud) test on sec%nb_point, initialisation of ztransp1,ztransp2,...
11   !!            3.4  ! 09/2011 (C Bricaud)
12   !!----------------------------------------------------------------------
13#if defined key_diadct
14   !!----------------------------------------------------------------------
15   !!   'key_diadct' :
16   !!----------------------------------------------------------------------
17   !!----------------------------------------------------------------------
18   !!   dia_dct      :  Compute the transport through a sec.
19   !!   dia_dct_init :  Read namelist.
20   !!   readsec      :  Read sections description and pathway
21   !!   removepoints :  Remove points which are common to 2 procs
22   !!   transport    :  Compute transport for each sections
23   !!   dia_dct_wri  :  Write tranports results in ascii files
24   !!   interp       :  Compute temperature/salinity/density at U-point or V-point
25   !!   
26   !!----------------------------------------------------------------------
27   USE oce             ! ocean dynamics and tracers
28   USE dom_oce         ! ocean space and time domain
29   USE phycst          ! physical constants
30   USE in_out_manager  ! I/O manager
31   USE daymod          ! calendar
32   USE dianam          ! build name of file
33   USE lib_mpp         ! distributed memory computing library
34#if defined key_si3
35   USE ice
36#endif
37   USE domvvl
38   USE timing          ! preformance summary
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   dia_dct      ! routine called by step.F90
44   PUBLIC   dia_dct_init ! routine called by opa.F90
45   PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90
46   PRIVATE  readsec
47   PRIVATE  removepoints
48   PRIVATE  transport
49   PRIVATE  dia_dct_wri
50
51   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
52
53   INTEGER :: nn_dct        ! Frequency of computation
54   INTEGER :: nn_dctwri     ! Frequency of output
55   INTEGER :: nn_secdebug   ! Number of the section to debug
56   
57   INTEGER, PARAMETER :: nb_class_max  = 10
58   INTEGER, PARAMETER :: nb_sec_max    = 150
59   INTEGER, PARAMETER :: nb_point_max  = 2000
60   INTEGER, PARAMETER :: nb_type_class = 10
61   INTEGER, PARAMETER :: nb_3d_vars    = 3 
62   INTEGER, PARAMETER :: nb_2d_vars    = 2 
63   INTEGER            :: nb_sec 
64
65   TYPE POINT_SECTION
66      INTEGER :: I,J
67   END TYPE POINT_SECTION
68
69   TYPE COORD_SECTION
70      REAL(wp) :: lon,lat
71   END TYPE COORD_SECTION
72
73   TYPE SECTION
74      CHARACTER(len=60)                            :: name              ! name of the sec
75      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
76                                                                       ! heat transports
77      LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
78      LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
79      TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
80      INTEGER                                      :: nb_class          ! number of boundaries for density classes
81      INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
82      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
83      REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
84                                                      zsigp           ,&! potential density classes    (99 if you don't want)
85                                                      zsal            ,&! salinity classes   (99 if you don't want)
86                                                      ztem            ,&! temperature classes(99 if you don't want)
87                                                      zlay              ! level classes      (99 if you don't want)
88      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
89      REAL(wp)                                         :: slopeSection  ! slope of the section
90      INTEGER                                          :: nb_point      ! number of points in the section
91      TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
92   END TYPE SECTION
93
94   TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
95 
96   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
97   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d 
98
99   !!----------------------------------------------------------------------
100   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
101   !! $Id$
102   !! Software governed by the CeCILL licence     (./LICENSE)
103   !!----------------------------------------------------------------------
104CONTAINS
105 
106  INTEGER FUNCTION diadct_alloc() 
107     !!----------------------------------------------------------------------
108     !!                   ***  FUNCTION diadct_alloc  ***
109     !!----------------------------------------------------------------------
110     INTEGER :: ierr(2) 
111     !!----------------------------------------------------------------------
112
113     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) ) 
114     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) ) 
115
116     diadct_alloc = MAXVAL( ierr ) 
117     IF( diadct_alloc /= 0 )   CALL ctl_warn('diadct_alloc: failed to allocate arrays') 
118 
119  END FUNCTION diadct_alloc 
120
121
122   SUBROUTINE dia_dct_init
123      !!---------------------------------------------------------------------
124      !!               ***  ROUTINE diadct  *** 
125      !!
126      !!  ** Purpose: Read the namelist parameters
127      !!              Open output files
128      !!
129      !!---------------------------------------------------------------------
130      INTEGER  ::   ios                 ! Local integer output status for namelist read
131      !!
132      NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
133      !!---------------------------------------------------------------------
134
135     REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections
136     READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
137901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp )
138
139     REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
140     READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
141902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp )
142     IF(lwm) WRITE ( numond, namdct )
143
144     IF( lwp ) THEN
145        WRITE(numout,*) " "
146        WRITE(numout,*) "diadct_init: compute transports through sections "
147        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
148        WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct
149        WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri
150
151        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
152                                            WRITE(numout,*)"       Debug section number: ", nn_secdebug 
153        ELSE IF ( nn_secdebug ==  0 )THEN ; WRITE(numout,*)"       No section to debug"
154        ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)"       Debug all sections"
155        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
156        ENDIF
157
158        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
159          &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
160
161     ENDIF
162
163     !Read section_ijglobal.diadct
164     CALL readsec
165
166     !open output file
167     IF( lwm ) THEN
168        CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
169        CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
170        CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
171     ENDIF
172
173     ! Initialise arrays to zero
174     transports_3d(:,:,:,:)=0.0 
175     transports_2d(:,:,:)  =0.0 
176     !
177  END SUBROUTINE dia_dct_init
178 
179 
180  SUBROUTINE dia_dct( kt )
181     !!---------------------------------------------------------------------
182     !!               ***  ROUTINE diadct  *** 
183     !!
184     !!  Purpose :: Compute section transports and write it in numdct files
185     !!   
186     !!  Method  :: All arrays initialised to zero in dct_init
187     !!             Each nn_dct time step call subroutine 'transports' for
188     !!               each section to sum the transports over each grid cell.
189     !!             Each nn_dctwri time step:
190     !!               Divide the arrays by the number of summations to gain
191     !!               an average value
192     !!               Call dia_dct_sum to sum relevant grid boxes to obtain
193     !!               totals for each class (density, depth, temp or sal)
194     !!               Call dia_dct_wri to write the transports into file
195     !!               Reinitialise all relevant arrays to zero
196     !!---------------------------------------------------------------------
197     INTEGER, INTENT(in) ::   kt
198     !
199     INTEGER ::   jsec              ! loop on sections
200     INTEGER ::   itotal            ! nb_sec_max*nb_type_class*nb_class_max
201     LOGICAL ::   lldebug =.FALSE.  ! debug a section 
202     INTEGER              , DIMENSION(1)    ::   ish     ! work array for mpp_sum
203     INTEGER              , DIMENSION(3)    ::   ish2    !   "
204     REAL(wp), ALLOCATABLE, DIMENSION(:)    ::   zwork   !   " 
205     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)::   zsum    !   "
206     !!---------------------------------------------------------------------   
207     !
208     IF( ln_timing )   CALL timing_start('dia_dct')
209
210     IF( lk_mpp )THEN
211        itotal = nb_sec_max*nb_type_class*nb_class_max
212        ALLOCATE( zwork(itotal) , zsum(nb_sec_max,nb_type_class,nb_class_max) )
213     ENDIF   
214 
215     ! Initialise arrays
216     zwork(:) = 0.0 
217     zsum(:,:,:) = 0.0
218
219     IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
220         WRITE(numout,*) " "
221         WRITE(numout,*) "diadct: compute transport"
222         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
223         WRITE(numout,*) "nb_sec = ",nb_sec
224     ENDIF
225
226 
227     ! Compute transport and write only at nn_dctwri
228     IF( MOD(kt,nn_dct)==0 ) THEN
229
230        DO jsec=1,nb_sec
231
232           !debug this section computing ?
233           lldebug=.FALSE.
234           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 ) lldebug=.TRUE. 
235
236           !Compute transport through section 
237           CALL transport(secs(jsec),lldebug,jsec) 
238
239        ENDDO
240             
241        IF( MOD(kt,nn_dctwri)==0 )THEN
242
243           IF( kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: average transports and write at kt = ",kt         
244 
245           !! divide arrays by nn_dctwri/nn_dct to obtain average
246           transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 
247           transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct) 
248 
249           ! Sum over each class
250           DO jsec=1,nb_sec 
251              CALL dia_dct_sum(secs(jsec),jsec) 
252           ENDDO 
253
254           !Sum on all procs
255           IF( lk_mpp )THEN
256              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
257              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
258              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
259              zwork(:)= RESHAPE(zsum(:,:,:), ish )
260              CALL mpp_sum(zwork, ish(1))
261              zsum(:,:,:)= RESHAPE(zwork,ish2)
262              DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
263           ENDIF
264
265           !Write the transport
266           DO jsec=1,nb_sec
267
268              IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
269           
270              !nullify transports values after writing
271              transports_3d(:,jsec,:,:)=0.
272              transports_2d(:,jsec,:  )=0.
273              secs(jsec)%transport(:,:)=0. 
274
275           ENDDO
276
277        ENDIF
278
279     ENDIF
280
281     IF( lk_mpp )THEN
282        itotal = nb_sec_max*nb_type_class*nb_class_max
283        DEALLOCATE( zwork , zsum  )
284     ENDIF   
285
286     IF( ln_timing )   CALL timing_stop('dia_dct')
287     !
288  END SUBROUTINE dia_dct
289
290
291  SUBROUTINE readsec 
292     !!---------------------------------------------------------------------
293     !!               ***  ROUTINE readsec  ***
294     !!
295     !!  ** Purpose:
296     !!            Read a binary file(section_ijglobal.diadct)
297     !!            generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
298     !!
299     !!
300     !!---------------------------------------------------------------------
301     INTEGER :: iptglo , iptloc                               ! Global and local number of points for a section
302     INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2  ! temporary  integer
303     INTEGER :: jsec, jpt                                     ! dummy loop indices
304     INTEGER, DIMENSION(2) :: icoord 
305     LOGICAL               :: llbon, lldebug   ! local logical
306     CHARACTER(len=160)    :: clname           ! filename
307     CHARACTER(len=200)    :: cltmp
308     CHARACTER(len=200)    :: clformat                          !automatic format
309     TYPE(POINT_SECTION),DIMENSION(nb_point_max)  ::coordtemp   !contains listpoints coordinates read in the file
310     INTEGER, DIMENSION(nb_point_max) :: directemp              !contains listpoints directions read in the files
311     !!-------------------------------------------------------------------------------------
312
313     !open input file
314     !---------------
315     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
316 
317     !---------------
318     !Read input file
319     !---------------
320     
321     DO jsec=1,nb_sec_max      !loop on the nb_sec sections
322
323        IF (  jsec==nn_secdebug .OR. nn_secdebug==-1  ) &
324           & WRITE(numout,*)'debuging for section number: ',jsec 
325
326        !initialization
327        !---------------
328        secs(jsec)%name=''
329        secs(jsec)%llstrpond    = .FALSE.  ; secs(jsec)%ll_ice_section = .FALSE.
330        secs(jsec)%ll_date_line = .FALSE.  ; secs(jsec)%nb_class       = 0
331        secs(jsec)%zsigi        = 99._wp   ; secs(jsec)%zsigp          = 99._wp
332        secs(jsec)%zsal         = 99._wp   ; secs(jsec)%ztem           = 99._wp
333        secs(jsec)%zlay         = 99._wp         
334        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0
335
336        !read section's number / name / computing choices / classes / slopeSection / points number
337        !-----------------------------------------------------------------------------------------
338        READ(numdct_in,iostat=iost)isec
339        IF (iost .NE. 0 )EXIT !end of file
340        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
341        IF( jsec .NE. isec )  CALL ctl_stop( cltmp )
342
343        IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )WRITE(numout,*)"isec ",isec 
344
345        READ(numdct_in)secs(jsec)%name
346        READ(numdct_in)secs(jsec)%llstrpond
347        READ(numdct_in)secs(jsec)%ll_ice_section
348        READ(numdct_in)secs(jsec)%ll_date_line
349        READ(numdct_in)secs(jsec)%coordSec
350        READ(numdct_in)secs(jsec)%nb_class
351        READ(numdct_in)secs(jsec)%zsigi
352        READ(numdct_in)secs(jsec)%zsigp
353        READ(numdct_in)secs(jsec)%zsal
354        READ(numdct_in)secs(jsec)%ztem
355        READ(numdct_in)secs(jsec)%zlay
356        READ(numdct_in)secs(jsec)%slopeSection
357        READ(numdct_in)iptglo
358
359        !debug
360        !-----
361
362        IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
363         
364            WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))' 
365
366            WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name)
367            WRITE(numout,*)       "      Compute heat and salt transport ? ",secs(jsec)%llstrpond
368            WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section
369            WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line
370            WRITE(numout,*)       "      Slope section :                   ",secs(jsec)%slopeSection
371            WRITE(numout,*)       "      Number of points in the section:  ",iptglo
372            WRITE(numout,*)       "      Number of classes                 ",secs(jsec)%nb_class
373            WRITE(numout,clformat)"      Insitu density classes :          ",secs(jsec)%zsigi
374            WRITE(numout,clformat)"      Potential density classes :       ",secs(jsec)%zsigp
375            WRITE(numout,clformat)"      Salinity classes :                ",secs(jsec)%zsal
376            WRITE(numout,clformat)"      Temperature classes :             ",secs(jsec)%ztem
377            WRITE(numout,clformat)"      Depth classes :                   ",secs(jsec)%zlay
378        ENDIF               
379
380        IF( iptglo /= 0 )THEN
381             
382           !read points'coordinates and directions
383           !--------------------------------------
384           coordtemp(:) = POINT_SECTION(0,0) !list of points read
385           directemp(:) = 0                  !value of directions of each points
386           DO jpt=1,iptglo
387              READ(numdct_in) i1, i2
388              coordtemp(jpt)%I = i1 
389              coordtemp(jpt)%J = i2
390           ENDDO
391           READ(numdct_in) directemp(1:iptglo)
392   
393           !debug
394           !-----
395           IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
396              WRITE(numout,*)"      List of points in global domain:"
397              DO jpt=1,iptglo
398                 WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt),directemp(jpt)
399              ENDDO                 
400           ENDIF
401 
402           !Now each proc selects only points that are in its domain:
403           !--------------------------------------------------------
404           iptloc = 0                    ! initialize number of points selected
405           DO jpt = 1, iptglo            ! loop on listpoint read in the file
406              !     
407              iiglo=coordtemp(jpt)%I          ! global coordinates of the point
408              ijglo=coordtemp(jpt)%J          !  "
409
410              IF( iiglo==jpiglo .AND. nimpp==1 )   iiglo = 2         !!gm BUG: Hard coded periodicity !
411
412              iiloc=iiglo-nimpp+1   ! local coordinates of the point
413              ijloc=ijglo-njmpp+1   !  "
414
415              !verify if the point is on the local domain:(1,nlei)*(1,nlej)
416              IF( iiloc >= 1 .AND. iiloc <= nlei .AND. &
417                  ijloc >= 1 .AND. ijloc <= nlej       )THEN
418                 iptloc = iptloc + 1                                                 ! count local points
419                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
420                 secs(jsec)%direction(iptloc) = directemp(jpt)                       ! store local direction
421              ENDIF
422              !
423           END DO
424     
425           secs(jsec)%nb_point=iptloc !store number of section's points
426
427           !debug
428           !-----
429           IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
430              WRITE(numout,*)"      List of points selected by the proc:"
431              DO jpt = 1,iptloc
432                 iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
433                 ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
434                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
435              ENDDO
436           ENDIF
437
438              IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
439              DO jpt = 1,iptloc
440                 iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
441                 ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
442              ENDDO
443              ENDIF
444
445           !remove redundant points between processors
446           !------------------------------------------
447           lldebug = .FALSE. ; IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) lldebug = .TRUE.
448           IF( iptloc .NE. 0 )THEN
449              CALL removepoints(secs(jsec),'I','top_list',lldebug)
450              CALL removepoints(secs(jsec),'I','bot_list',lldebug)
451              CALL removepoints(secs(jsec),'J','top_list',lldebug)
452              CALL removepoints(secs(jsec),'J','bot_list',lldebug)
453           ENDIF
454           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
455              DO jpt = 1,secs(jsec)%nb_point
456                 iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
457                 ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
458              ENDDO
459           ENDIF
460
461           !debug
462           !-----
463           IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
464              WRITE(numout,*)"      List of points after removepoints:"
465              iptloc = secs(jsec)%nb_point
466              DO jpt = 1,iptloc
467                 iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
468                 ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
469                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
470                 CALL FLUSH(numout)
471              ENDDO
472           ENDIF
473
474        ELSE  ! iptglo = 0
475           IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )&
476              WRITE(numout,*)'   No points for this section.'
477        ENDIF
478
479     ENDDO !end of the loop on jsec
480 
481     nb_sec = jsec-1   !number of section read in the file
482     !
483  END SUBROUTINE readsec
484
485
486  SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
487     !!---------------------------------------------------------------------------
488     !!             *** function removepoints
489     !!
490     !!   ** Purpose :: Remove points which are common to 2 procs
491     !!
492     !----------------------------------------------------------------------------
493     TYPE(SECTION),INTENT(INOUT) :: sec
494     CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
495     CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
496     LOGICAL,INTENT(IN)          :: ld_debug                     
497     !
498     INTEGER :: iextr         ,& !extremity of listpoint that we verify
499                iind          ,& !coord     of listpoint that we verify
500                itest         ,& !indice value of the side of the domain
501                                 !where points could be redundant
502                isgn          ,& ! isgn= 1 : scan listpoint from start to end
503                                 ! isgn=-1 : scan listpoint from end to start
504                istart,iend      !first and last points selected in listpoint
505     INTEGER :: jpoint           !loop on list points
506     INTEGER, DIMENSION(nb_point_max)   :: idirec !contains temporary sec%direction
507     INTEGER, DIMENSION(2,nb_point_max) :: icoord !contains temporary sec%listpoint
508     !----------------------------------------------------------------------------
509     !
510     IF( ld_debug )WRITE(numout,*)'      -------------------------'
511     IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
512
513     !iextr=extremity of list_point that we verify
514     IF      ( cdextr=='bot_list' )THEN ; iextr=1            ; isgn=1
515     ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
516     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdextr")
517     ENDIF
518 
519     !which coordinate shall we verify ?
520     IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1
521     ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2
522     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind") 
523     ENDIF
524
525     IF( ld_debug )THEN
526        WRITE(numout,*)'      case: coord/list extr/domain side'
527        WRITE(numout,*)'      ', cdind,' ',cdextr,' ',itest
528        WRITE(numout,*)'      Actual number of points: ',sec%nb_point
529     ENDIF
530
531     icoord(1,1:nb_point_max) = sec%listPoint%I
532     icoord(2,1:nb_point_max) = sec%listPoint%J
533     idirec                   = sec%direction
534     sec%listPoint            = POINT_SECTION(0,0)
535     sec%direction            = 0
536
537     jpoint=iextr+isgn
538     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
539         IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
540         ELSE                                                                               ; EXIT
541         ENDIF
542     ENDDO 
543
544     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
545     ELSE                        ; istart=1        ; iend=jpoint+1
546     ENDIF
547
548     sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
549     sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
550     sec%direction(1:1+iend-istart)   = idirec(istart:iend)
551     sec%nb_point                     = iend-istart+1
552     
553     IF( ld_debug )THEN
554        WRITE(numout,*)'      Number of points after removepoints :',sec%nb_point
555        WRITE(numout,*)'      sec%direction after removepoints :',sec%direction(1:sec%nb_point)
556     ENDIF
557      !
558   END SUBROUTINE removepoints
559
560
561   SUBROUTINE transport(sec,ld_debug,jsec)
562     !!-------------------------------------------------------------------------------------------
563     !!                     ***  ROUTINE transport  ***
564     !!
565     !!  Purpose ::  Compute the transport for each point in a section
566     !!
567     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
568     !!              Be aware :           
569     !!              One section is a sum of segments
570     !!              One segment is defined by 2 consecutive points in sec%listPoint
571     !!              All points of sec%listPoint are positioned on the F-point of the cell
572     !!
573     !!              There are two loops:                 
574     !!              loop on the segment between 2 nodes
575     !!              loop on the level jk !!
576     !!
577     !!  Output  ::  Arrays containing the volume,density,heat,salt transports for each i
578     !!              point in a section, summed over each nn_dct.
579     !!
580     !!-------------------------------------------------------------------------------------------
581     TYPE(SECTION),INTENT(INOUT) :: sec
582     LOGICAL      ,INTENT(IN)    :: ld_debug
583     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
584     !
585     INTEGER ::   jk, jseg, jclass,jl, isgnu, isgnv    ! loop on level/segment/classes/ice categories
586     REAL(wp)::   zumid, zvmid, zumid_ice, zvmid_ice   ! U/V ocean & ice velocity on a cell segment
587     REAL(wp)::   zTnorm                               ! transport of velocity through one cell's sides
588     REAL(wp)::   ztn, zsn, zrhoi, zrhop, zssh, zdep   ! temperature/salinity/potential density/ssh/depth at u/v point
589     TYPE(POINT_SECTION) ::   k
590      !!--------------------------------------------------------
591      !
592      IF( ld_debug )WRITE(numout,*)'      Compute transport'
593
594      !---------------------------!
595      !  COMPUTE TRANSPORT        !
596      !---------------------------!
597      IF(sec%nb_point .NE. 0)THEN   
598
599         !----------------------------------------------------------------------------------------------------
600         !Compute sign for velocities:
601         !
602         !convention:
603         !   non horizontal section: direction + is toward left hand of section
604         !       horizontal section: direction + is toward north of section
605         !
606         !
607         !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
608         !       ----------------      -----------------     ---------------             --------------
609         !
610         !   isgnv=1         direction +     
611         !  ______         _____             ______                                                   
612         !        |           //|            |                  |                         direction +   
613         !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
614         !        |_______  //         ______|    \\            | ---\                        |
615         !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
616         !               |             |          __\\|         |                   
617         !               |             |     direction +        |                      isgnv=1                                 
618         !                                                     
619         !----------------------------------------------------------------------------------------------------
620         isgnu = 1
621         IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
622         ELSE                                ; isgnv =  1
623         ENDIF
624         IF( sec%slopeSection .GE. 9999. )     isgnv =  1
625
626         IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
627
628         !--------------------------------------!
629         ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
630         !--------------------------------------!
631         DO jseg=1,MAX(sec%nb_point-1,0)
632             
633            !-------------------------------------------------------------------------------------------
634            ! Select the appropriate coordinate for computing the velocity of the segment
635            !
636            !                      CASE(0)                                    Case (2)
637            !                      -------                                    --------
638            !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
639            !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
640            !                                                                            |
641            !                                                                            |
642            !                                                                            |
643            !                      Case (3)                                            U(i,j)
644            !                      --------                                              |
645            !                                                                            |
646            !  listPoint(jseg+1) F(i,j+1)                                                |
647            !                        |                                                   |
648            !                        |                                                   |
649            !                        |                                 listPoint(jseg+1) F(i,j-1)
650            !                        |                                           
651            !                        |                                           
652            !                     U(i,j+1)                                           
653            !                        |                                       Case(1)     
654            !                        |                                       ------     
655            !                        |                                           
656            !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
657            !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
658            ! listPoint(jseg)     F(i,j)
659            !
660            !-------------------------------------------------------------------------------------------
661
662            SELECT CASE( sec%direction(jseg) )
663            CASE(0)   ;    k = sec%listPoint(jseg)
664            CASE(1)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
665            CASE(2)   ;    k = sec%listPoint(jseg)
666            CASE(3)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
667            END SELECT
668
669            !---------------------------|
670            !     LOOP ON THE LEVEL     |
671            !---------------------------|
672            DO jk = 1, mbkt(k%I,k%J)            !Sum of the transport on the vertical
673            !           ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
674            SELECT CASE( sec%direction(jseg) )
675               CASE(0,1) 
676                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
677                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
678                  zrhop = interp(k%I,k%J,jk,'V',rhop) 
679                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0) 
680                  zssh  =  0.5*( ssh(k%I,k%J,Nnn) + ssh(k%I,k%J+1,Nnn)    ) * vmask(k%I,k%J,1) 
681               CASE(2,3) 
682                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
683                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
684                  zrhop = interp(k%I,k%J,jk,'U',rhop) 
685                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0) 
686                  zssh  =  0.5*( ssh(k%I,k%J,Nnn) + ssh(k%I+1,k%J,Nnn)    ) * umask(k%I,k%J,1) 
687               END SELECT 
688               !
689               zdep= gdept_n(k%I,k%J,jk) 
690 
691               SELECT CASE( sec%direction(jseg) )                !compute velocity with the correct direction
692               CASE(0,1)   
693                  zumid=0._wp
694                  zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
695               CASE(2,3) 
696                  zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
697                  zvmid=0._wp
698               END SELECT 
699 
700               !zTnorm=transport through one cell;
701               !velocity* cell's length * cell's thickness
702               zTnorm = zumid*e2u(k%I,k%J) * e3u_n(k%I,k%J,jk)     & 
703                  &   + zvmid*e1v(k%I,k%J) * e3v_n(k%I,k%J,jk) 
704
705!!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!!
706!               IF( ln_linssh ) THEN              !add transport due to free surface
707!                  IF( jk==1 ) THEN
708!                     zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zssh * umask(k%I,k%J,jk)   &
709!                        &            + zvmid* e1v(k%I,k%J) * zssh * vmask(k%I,k%J,jk)
710!                  ENDIF
711!               ENDIF
712!!gm end
713              !COMPUTE TRANSPORT 
714 
715              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
716 
717              IF( sec%llstrpond ) THEN
718                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp
719                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001
720              ENDIF
721   
722           END DO !end of loop on the level
723
724#if defined key_si3
725
726           !ICE CASE   
727           !------------
728           IF( sec%ll_ice_section )THEN
729              SELECT CASE (sec%direction(jseg))
730              CASE(0)
731                 zumid_ice = 0
732                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
733              CASE(1)
734                 zumid_ice = 0
735                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
736              CASE(2)
737                 zvmid_ice = 0
738                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
739              CASE(3)
740                 zvmid_ice = 0
741                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
742              END SELECT
743   
744              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
745
746#if defined key_si3
747              DO jl=1,jpl
748                 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*       &
749                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) *  &
750                                  ( h_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  &
751                                    h_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
752                                   
753                 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
754                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
755              END DO
756#endif
757   
758           ENDIF !end of ice case
759#endif
760 
761        END DO !end of loop on the segment
762
763     ENDIF !end of sec%nb_point =0 case
764     !
765  END SUBROUTINE transport
766
767
768  SUBROUTINE dia_dct_sum(sec,jsec) 
769     !!-------------------------------------------------------------
770     !! Purpose: Average the transport over nn_dctwri time steps 
771     !! and sum over the density/salinity/temperature/depth classes
772     !!
773     !! Method:   Sum over relevant grid cells to obtain values 
774     !!           for each class
775     !!              There are several loops:                 
776     !!              loop on the segment between 2 nodes
777     !!              loop on the level jk
778     !!              loop on the density/temperature/salinity/level classes
779     !!              test on the density/temperature/salinity/level
780     !!
781     !!  Note:    Transport through a given section is equal to the sum of transports
782     !!           computed on each proc.
783     !!           On each proc,transport is equal to the sum of transport computed through
784     !!           segments linking each point of sec%listPoint  with the next one.   
785     !!
786     !!-------------------------------------------------------------
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, zssh, 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(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
849                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
850                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
851                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0) 
852
853              CASE(2,3) 
854                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
855                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
856                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
857                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0) 
858                 zssh  =  0.5*(  ssh(k%I,k%J,Nnn) + ssh(k%I+1,k%J,Nnn)  ) * umask(k%I,k%J,1) 
859              END SELECT
860 
861              zdep= gdept_n(k%I,k%J,jk) 
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(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)                          :: ki, kj, kk   ! coordinate of point
1166  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
1167  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk )
1168  REAL(wp)                                     :: interp       ! interpolated variable
1169
1170  !*local declations
1171  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
1172  REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu            ! local real
1173  REAL(wp):: zet1, zet2                                        ! weight for interpolation
1174  REAL(wp):: zdep1,zdep2                                       ! differences of depth
1175  REAL(wp):: zmsk                                              ! mask value
1176  !!----------------------------------------------------------------------
1177
1178  IF( cd_point=='U' )THEN
1179     ii1 = ki    ; ij1 = kj 
1180     ii2 = ki+1  ; ij2 = kj 
1181
1182     zet1=e1t(ii1,ij1)
1183     zet2=e1t(ii2,ij2)
1184     zmsk=umask(ii1,ij1,kk)
1185 
1186
1187  ELSE ! cd_point=='V'
1188     ii1 = ki    ; ij1 = kj 
1189     ii2 = ki    ; ij2 = kj+1 
1190
1191     zet1=e2t(ii1,ij1)
1192     zet2=e2t(ii2,ij2)
1193     zmsk=vmask(ii1,ij1,kk)
1194
1195  ENDIF
1196
1197  IF( ln_sco )THEN   ! s-coordinate case
1198
1199     zdepu = ( gdept_n(ii1,ij1,kk) +  gdept_n(ii2,ij2,kk) ) * 0.5_wp 
1200     zdep1 = gdept_n(ii1,ij1,kk) - zdepu
1201     zdep2 = gdept_n(ii2,ij2,kk) - zdepu
1202
1203     ! weights
1204     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1205     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1206 
1207     ! result
1208     interp = zmsk * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )   
1209
1210
1211  ELSE       ! full step or partial step case
1212
1213     ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk) 
1214     zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
1215     zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
1216
1217     IF(kk .NE. 1)THEN
1218
1219        IF( ze3t >= 0. )THEN 
1220           ! zbis
1221           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
1222           ! result
1223            interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
1224        ELSE
1225           ! zbis
1226           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1227           ! result
1228           interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1229        ENDIF   
1230
1231     ELSE
1232        interp = zmsk * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1233     ENDIF
1234
1235  ENDIF
1236      !
1237   END FUNCTION interp
1238
1239#else
1240   !!----------------------------------------------------------------------
1241   !!   Default option :                                       Dummy module
1242   !!----------------------------------------------------------------------
1243   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
1244   PUBLIC 
1245   !! $Id$
1246CONTAINS
1247
1248   SUBROUTINE dia_dct_init          ! Dummy routine
1249      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
1250   END SUBROUTINE dia_dct_init
1251
1252   SUBROUTINE dia_dct( kt )         ! Dummy routine
1253      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
1254      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
1255   END SUBROUTINE dia_dct
1256#endif
1257
1258   !!======================================================================
1259END MODULE diadct
Note: See TracBrowser for help on using the repository browser.