source: branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 2864

Last change on this file since 2864 was 2864, checked in by cbricaud, 9 years ago

add dia_dct_init Dummy routine


  • Property svn:executable set to *
File size: 47.0 KB
Line 
1MODULE diadct
2  !!=====================================================================
3  !!                       ***  MODULE  diadct  ***
4  !! Ocean diagnostics: Compute the transport trough a sec.
5  !!===============================================================
6  !! History :
7  !!
8  !!         original  : 02/99 (Y Drillet)
9  !!         addition  : 10/01 (Y Drillet, R Bourdalle Badie)
10  !!                   : 10/05 (M Laborie) F90
11  !!         addition  : 04/07 (G Garric) Ice sections
12  !!         bugfix    : 04/07 (C Bricaud) test on sec%nb_point
13  !!                                      initialisation of ztransp1,ztransp2,...
14  !!         nemo_v_3_4: 09/2011 (C Bricaud)
15  !!
16  !!
17  !!----------------------------------------------------------------------
18#if defined key_diadct
19  !!----------------------------------------------------------------------
20  !!   'key_diadct' :
21  !!----------------------------------------------------------------------
22  !!----------------------------------------------------------------------
23  !!   dia_dct      :  compute the transport trough a sec.
24  !!   dia_dct_init :  read namelist.
25  !!   readsec      :  read sections description and pathway
26  !!   removepoints :  remove points which are common to 2 procs
27  !!   transport    :  Compute transport for each sections
28  !!   dia_dct_wri  :  write tranports results in ascii files
29  !!   interp       :  compute Temperature/Salinity/density on U-point or V-point
30  !!   
31  !!----------------------------------------------------------------------
32  !! * Modules used
33  USE oce             ! ocean dynamics and tracers
34  USE dom_oce         ! ocean space and time domain
35  USE phycst          ! physical constants
36  USE in_out_manager  ! I/O manager
37  USE daymod          ! calendar
38  USE dianam          ! build name of file
39  USE lib_mpp         ! distributed memory computing library
40#if defined key_ice_lim2 || defined key_ice_lim3
41  USE ice
42#endif
43
44  IMPLICIT NONE
45  PRIVATE
46
47  !! * Routine accessibility
48  PUBLIC   dia_dct     ! routine called by step.F90
49  PUBLIC   dia_dct_init! routine called by opa.F90
50  PRIVATE  readsec
51  PRIVATE  removepoints
52  PRIVATE  transport
53  PRIVATE  dia_dct_wri
54
55#include "domzgr_substitute.h90"
56
57  !! * Shared module variables
58  LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
59
60  !! * Module variables
61  INTEGER :: nn_dct      = 1     ! Frequency of computation
62  INTEGER :: nn_dctwri   = 1     ! Frequency of output
63  INTEGER :: nn_secdebug = 0     ! Number of the section to debug
64   
65  INTEGER, PARAMETER :: nb_class_max  = 10
66  INTEGER, PARAMETER :: nb_sec_max    = 150
67  INTEGER, PARAMETER :: nb_point_max  = 2000
68  INTEGER, PARAMETER :: nb_type_class = 14
69  INTEGER            :: nb_sec 
70 
71  TYPE POINT_SECTION
72     INTEGER :: I,J
73  END TYPE POINT_SECTION
74
75  TYPE COORD_SECTION
76     REAL(wp) :: lon,lat
77  END TYPE COORD_SECTION
78
79  TYPE SECTION
80     CHARACTER(len=60)                            :: name           ! name of the sec
81     LOGICAL                                      :: llstrpond      ! true if you want the computation of salinity and
82                                                                    ! temperature balanced by the transport
83     LOGICAL                                      :: ll_ice_section ! icesurf and icevol computation
84     LOGICAL                                      :: ll_date_line   ! = T if the section crosses the date-line
85     TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec       ! longitude and latitude of the extremities of the sec
86     INTEGER                                      :: nb_class       ! number of boundaries for density classes
87     INTEGER, DIMENSION(nb_point_max)             :: direction      ! vector direction of the point in the section
88     CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname      ! caracteristics of the class
89     REAL(wp), DIMENSION(nb_class_max)            :: zsigi        ,&! in-situ   density classes    (99 if you don't want)
90                                                     zsigp        ,&! potential density classes    (99 if you don't want)
91                                                     zsal         ,&! salinity classes   (99 if you don't want)
92                                                     ztem         ,&! temperature classes(99 if you don't want)
93                                                     zlay           ! level classes      (99 if you don't want)
94     REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
95     REAL(wp)                                         :: slopeSection  ! coeff directeur de la section
96     INTEGER                                          :: nb_point              ! nb de points de la section
97     TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint             ! list of point in the section
98  END TYPE SECTION
99
100  TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
101 
102 
103CONTAINS
104
105  SUBROUTINE dia_dct_init
106     !!---------------------------------------------------------------------
107     !!               ***  ROUTINE diadct  *** 
108     !!
109     !!  ** Purpose: Read the namelist parametres
110     !!              Open output files
111     !!
112     !!---------------------------------------------------------------------
113     NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
114
115     !read namelist
116     REWIND( numnam )
117     READ  ( numnam, namdct )
118
119     IF( lwp ) THEN
120        WRITE(numout,*) " "
121        WRITE(numout,*) "diadct_init: compute transports through sections "
122        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
123        WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct
124        WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri
125
126        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
127                                            WRITE(numout,*)"       Debug section number: ", nn_secdebug 
128        ELSE IF ( nn_secdebug ==  0 )THEN ; WRITE(numout,*)"       No section to debug"
129        ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)"       Debug all sections"
130        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
131        ENDIF
132
133        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
134          &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
135
136        !open output file
137        CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
138        CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
139        CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
140
141     ENDIF
142
143     !Read section_ijglobal.diadct
144     CALL readsec
145
146  END SUBROUTINE dia_dct_init
147 
148 
149  SUBROUTINE dia_dct(kt)
150     !!---------------------------------------------------------------------
151     !!               ***  ROUTINE diadct  *** 
152     !!
153     !!  ** Purpose: Compute sections tranport and write it in numdct file
154     !!---------------------------------------------------------------------
155     !! * Arguments
156     INTEGER,INTENT(IN)        ::kt
157
158     !! * Local variables
159     INTEGER             :: jsec,            &!loop on sections
160                            iost              !error for opening fileout
161     LOGICAL             :: lldebug =.FALSE.  !debug a section 
162     CHARACTER(len=160)  :: clfileout         !fileout name
163     !!---------------------------------------------------------------------   
164
165     IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
166         WRITE(numout,*) " "
167         WRITE(numout,*) "diadct: compute transport"
168         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
169         WRITE(numout,*) "nb_sec = ",nb_sec
170     ENDIF
171
172 
173     ! Compute transport and write only at nn_dctwri
174     IF( MOD(kt,nn_dct)==0 ) THEN
175
176        DO jsec=1,nb_sec
177
178           !debug this section computing ?
179           lldebug=.FALSE.
180!           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE.
181           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 ) lldebug=.TRUE. 
182
183           !Compute transport through section 
184           CALL transport(secs(jsec),lldebug) 
185             
186           IF( MOD(kt,nn_dctwri)==0 )THEN
187
188              IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: write at kt = ",kt         
189 
190              !Write the transport
191              IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec))
192           
193              !nullify transports values after writing
194              secs(jsec)%transport(:,:)=0. 
195
196           ENDIF
197        ENDDO
198     ENDIF
199
200  END SUBROUTINE dia_dct
201
202  SUBROUTINE readsec 
203     !!---------------------------------------------------------------------
204     !!               ***  ROUTINE readsec  ***
205     !!
206     !!  ** Purpose:
207     !!            Read a binary file(section_ijglobal.diadct)
208     !!            generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
209     !!
210     !!
211     !!---------------------------------------------------------------------
212     !! * Local variables
213     INTEGER :: iptglo , iptloc                               ! Global and local number of points for a section
214     INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2  ! temporary  integer
215     INTEGER :: jsec, jpt                                     ! dummy loop indices
216
217     INTEGER, DIMENSION(2) :: icoord 
218     CHARACTER(len=160)  :: clname                            !filename
219     CHARACTER(len=200)  :: cltmp
220     TYPE(POINT_SECTION),DIMENSION(nb_point_max)  ::coordtemp !contains listpoints coordinates
221                                                              !read in the file
222     INTEGER,DIMENSION(nb_point_max)  ::directemp             !contains listpoints directions
223                                                              !read in the files
224     LOGICAL :: llbon                                       ,&!local logical
225                lldebug                                       !debug the section
226     !!-------------------------------------------------------------------------------------
227
228     !open input file
229     !---------------
230     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
231 
232     !---------------
233     !Read input file
234     !---------------
235     
236     DO jsec=1,nb_sec_max      !loop on the nb_sec sections
237
238        IF ( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) ) &
239           & WRITE(numout,*)'debuging for section number: ',jsec 
240
241        !initialization
242        !---------------
243        secs(jsec)%name=''
244        secs(jsec)%llstrpond    = .FALSE.  ; secs(jsec)%ll_ice_section = .FALSE.
245        secs(jsec)%ll_date_line = .FALSE.  ; secs(jsec)%nb_class       = 0
246        secs(jsec)%zsigi        = 99._wp   ; secs(jsec)%zsigp          = 99._wp
247        secs(jsec)%zsal         = 99._wp   ; secs(jsec)%ztem           = 99._wp
248        secs(jsec)%zlay         = 99._wp         
249        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0
250
251        !read section's number / name / computing choices / classes / slopeSection / points number
252        !-----------------------------------------------------------------------------------------
253        READ(numdct_in,iostat=iost)isec
254        IF (iost .NE. 0 )EXIT !end of file
255        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' et jsec= ',jsec
256        IF( jsec .NE. isec )  CALL ctl_stop( cltmp )
257
258        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec 
259
260        READ(numdct_in)secs(jsec)%name
261        READ(numdct_in)secs(jsec)%llstrpond
262        READ(numdct_in)secs(jsec)%ll_ice_section
263        READ(numdct_in)secs(jsec)%ll_date_line
264        READ(numdct_in)secs(jsec)%coordSec
265        READ(numdct_in)secs(jsec)%nb_class
266        READ(numdct_in)secs(jsec)%zsigi
267        READ(numdct_in)secs(jsec)%zsigp
268        READ(numdct_in)secs(jsec)%zsal
269        READ(numdct_in)secs(jsec)%ztem
270        READ(numdct_in)secs(jsec)%zlay
271        READ(numdct_in)secs(jsec)%slopeSection
272        READ(numdct_in)iptglo
273
274        !debug
275        !-----
276        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
277            WRITE(numout,*)                  "   Section name :                       ",TRIM(secs(jsec)%name)
278            WRITE(numout,*)                  "      Compute heat and salt transport ? ",secs(jsec)%llstrpond
279            WRITE(numout,*)                  "      Compute ice transport ?           ",secs(jsec)%ll_ice_section
280            WRITE(numout,*)                  "      Section crosses date-line ?       ",secs(jsec)%ll_date_line
281            WRITE(numout,*)                  "      Slope section :                   ",secs(jsec)%slopeSection
282            WRITE(numout,*)                  "      Number of points in the section:  ",iptglo
283            WRITE(numout,*)                  "      Number of classes                 ",secs(jsec)%nb_class
284            WRITE(numout,'(A40,10(f8.3,1X))')"      Insitu density classes :          ",secs(jsec)%zsigi
285            WRITE(numout,'(A40,10(f8.3,1X))')"      Potential density classes :       ",secs(jsec)%zsigp
286            WRITE(numout,'(A40,10(f8.3,1X))')"      Salinity classes :                ",secs(jsec)%zsal
287            WRITE(numout,'(A40,10(f8.3,1X))')"      Temperature classes :             ",secs(jsec)%ztem
288            WRITE(numout,'(A40,10(f8.3,1X))')"      Depth classes :                   ",secs(jsec)%zlay
289        ENDIF               
290
291        IF( iptglo .NE. 0 )THEN
292             
293           !read points'coordinates and directions
294           !--------------------------------------
295           coordtemp(:) = POINT_SECTION(0,0) !list of points read
296           directemp(:) = 0                  !value of directions of each points
297           DO jpt=1,iptglo
298              READ(numdct_in)i1,i2
299              coordtemp(jpt)%I = i1 
300              coordtemp(jpt)%J = i2
301           ENDDO
302           READ(numdct_in)directemp(1:iptglo)
303   
304           !debug
305           !-----
306           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
307              WRITE(numout,*)"      List of points in global domain:"
308              DO jpt=1,iptglo
309                 WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt)
310              ENDDO                 
311           ENDIF
312           !?IF( jsec==nn_secdebug .OR. nn_secdebug==-1  )THEN
313           !  DO jpt=1,iptglo
314           !      WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt)
315           !   ENDDO                 
316           !ENDIF
317 
318           !Now each proc selects only points that are in its domain:
319           !--------------------------------------------------------
320           iptloc = 0                    !initialize number of points selected
321           DO jpt=1,iptglo               !loop on listpoint read in the file
322                   
323              iiglo=coordtemp(jpt)%I          ! global coordinates of the point
324              ijglo=coordtemp(jpt)%J          !  "
325
326              IF( iiglo==jpidta .AND. nimpp==1 ) iiglo = 2
327
328              iiloc=iiglo-jpizoom+1-nimpp+1   ! local coordinates of the point
329              ijloc=ijglo-jpjzoom+1-njmpp+1   !  "
330
331              !verify if the point is on the local domain:(1,nlei)*(1,nlej)
332              IF( iiloc .GE. 1 .AND. iiloc .LE. nlei .AND. &
333                  ijloc .GE. 1 .AND. ijloc .LE. nlej       )THEN
334                 iptloc = iptloc + 1                                                 ! count local points
335                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
336                 secs(jsec)%direction(iptloc) = directemp(jpt)               ! store local direction
337                 !WRITE(narea+200,*)'         # I J : ',iiloc,ijloc
338              ENDIF
339
340           ENDDO
341     
342           secs(jsec)%nb_point=iptloc !store number of section's points
343
344           !debug
345           !-----
346           !IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
347           IF(  ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
348              WRITE(numout,*)"      List of points selected by the proc:"
349              DO jpt = 1,iptloc
350                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
351                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
352                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
353                 !WRITE(narea+200,*)'         # I J : ',iiglo,ijglo
354              ENDDO
355           ENDIF
356
357           !remove redundant points between processors
358           !look NEMO documentation for more explanations
359!           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
360           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1)  ) lldebug = .TRUE.
361           IF( iptloc .NE. 0 )THEN
362              CALL removepoints(secs(jsec),'I','top_list',lldebug)
363              CALL removepoints(secs(jsec),'I','bot_list',lldebug)
364              CALL removepoints(secs(jsec),'J','top_list',lldebug)
365              CALL removepoints(secs(jsec),'J','bot_list',lldebug)
366           ENDIF
367
368           !debug
369           !-----
370           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
371              WRITE(numout,*)"      List of points after removepoints:"
372              DO jpt = 1,iptloc
373                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
374                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
375                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
376              ENDDO
377           ENDIF
378
379        ELSE  ! iptglo = 0
380           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )&
381              WRITE(numout,*)'   Not points for this section.'
382        ENDIF
383
384     ENDDO !end of the loop on jsec
385 
386     nb_sec = jsec-1   !number of section read in the file
387
388  END SUBROUTINE readsec
389
390  SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
391     !!---------------------------------------------------------------------------
392     !!             *** function removepoints
393     !!
394     !!   ** Purpose ::
395     !!              remove points which are common to 2 procs
396     !!
397     !!
398     !----------------------------------------------------------------------------
399     !! * arguments
400     TYPE(SECTION),INTENT(INOUT) :: sec
401     CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
402     CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
403     LOGICAL,INTENT(IN)          :: ld_debug                     
404
405     !! * Local variables
406     INTEGER :: iextr         ,& !extremity of listpoint that we verify
407                iind          ,& !coord     of listpoint that we verify
408                itest         ,& !indice value of the side of the domain
409                                 !where points could be redundant
410                isgn          ,& !way of course in listpoint
411                istart,iend      !first and last points selected in listpoint
412     INTEGER :: jpoint   =0      !loop on list points
413     INTEGER,DIMENSION(nb_point_max)   :: idirec !contains temporare sec%direction
414     INTEGER,DIMENSION(2,nb_point_max) :: icoord !contains temporare sec%listpoint
415     !----------------------------------------------------------------------------
416     IF( ld_debug )WRITE(numout,*)'      -------------------------'
417     IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
418
419     !iextr=extremity of list_point that we verify
420     IF      ( cdextr=='bot_list' )THEN ; iextr=1            ; isgn=1
421     ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
422     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdextr")
423     ENDIF
424 
425     !which coordinate shall we verify ?
426     IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1
427     ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2
428     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind") 
429     ENDIF
430
431     IF( ld_debug )THEN
432        WRITE(numout,*)'      case: coord/list extr/domain side'
433        WRITE(numout,*)'      ', cdind,' ',cdextr,' ',itest
434        WRITE(numout,*)'      Actual number of points: ',sec%nb_point
435     ENDIF
436
437     icoord(1,1:nb_point_max) = sec%listPoint%I
438     icoord(2,1:nb_point_max) = sec%listPoint%J
439     idirec                   = sec%direction
440     sec%listPoint            = POINT_SECTION(0,0)
441     sec%direction            = 0
442
443
444     jpoint=iextr+isgn
445     DO WHILE ( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point .AND. &
446             icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest)
447             jpoint=jpoint+isgn
448     ENDDO
449     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
450     ELSE                        ; istart=1        ; iend=jpoint+1
451     ENDIF
452     sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
453     sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
454     sec%direction(1:1+iend-istart)   = idirec(istart:iend)
455     sec%nb_point                     = iend-istart+1
456     
457     IF( ld_debug )THEN
458        WRITE(numout,*)'      Number of points after removepoints :',sec%nb_point
459        WRITE(numout,*)'      sec%direction after removepoints :',sec%direction(1:sec%nb_point)
460     ENDIF
461
462  END SUBROUTINE removepoints
463
464  SUBROUTINE transport(sec,ld_debug)
465     !!---------------------------------------------------------------------
466     !!                     ***  ROUTINE transport  ***
467     !!
468     !!  ** Purpose : Compute the transport trough a sec
469     !!
470     !!  ** Method  :Transport through a given section is equal to the sum of transports
471     !!              computed on each proc.
472     !!              On each proc,transport is equal to the sum of transport computed through
473     !!               segments linking each points of sec%listPoint  with the next one.           
474     !!              There are several loops:                 
475     !!              loop on the density/temperature/salinity/level classes
476     !!              loop on the segment between 2 nodes
477     !!              loop on the level jk
478     !!              test on the density/temperature/salinity/level
479     !!
480     !! ** Output: sec%transport: transport in the 2 direction and temperature,
481     !!                       salinity, density meaned by the transport
482     !!
483     !!
484     !!-------------------------------------------------------------------
485     !! * Arguments
486     TYPE(SECTION),INTENT(INOUT) :: sec
487     LOGICAL      ,INTENT(IN)    :: ld_debug
488   
489     !! * Local variables
490     INTEGER             :: jk,jseg,jclass,   &!loop on level/segment/classes
491                            iptegrid,         &!=0 if section is constant along i/j
492                            isgnu  , isgnv     !
493     INTEGER :: ii, ij ! local integer
494     REAL(wp):: zumid        , zvmid        ,&!U/V velocity on a cell segment
495                zumid_ice    , zvmid_ice    ,&!U/V ice velocity
496                zTnorm                      ,&!transport of velocity through one cell's sides
497                ztransp1     , ztransp2     ,&!total        transport in sens 1 and 2
498                ztemp1       , ztemp2       ,&!temperature  transport     "
499                zrhoi1       , zrhoi2       ,&!mass         transport     "
500                zrhop1       , zrhop2       ,&!mass         transport     "
501                zsal1        , zsal2        ,&!salinity     transport     "
502                zice_vol_pos , zice_vol_neg ,&!volumic  ice transport     "
503                zice_surf_pos, zice_surf_neg  !surfacic ice transport     "
504     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
505     REAL(wp) :: aj,i0,j0,i1,j1,i,j
506
507     TYPE(POINT_SECTION) :: k
508     INTEGER,DIMENSION(1):: ish
509     INTEGER,DIMENSION(2):: ish2
510     REAL(wp),DIMENSION(nb_type_class,nb_class_max)::zsum
511     REAL(wp),DIMENSION(nb_type_class*nb_class_max)::zwork  ! temporary vector for mpp_sum
512     !!--------------------------------------------------------
513
514     IF( ld_debug )WRITE(numout,*)'      Compute transport'
515
516     !----------------!
517     ! INITIALIZATION !
518     !----------------!
519     zsum    = 0._wp
520     zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp
521     zice_vol_pos  = 0._wp ; zice_vol_neg  = 0._wp
522
523     !compute iptegrid: iptegrid=0 if section is constant along i or j
524     iptegrid = 0
525     jseg     = 1
526     DO WHILE( (jseg .LT. MAX(sec%nb_point-1,0)) .AND. iptegrid .NE. 1 )
527        IF(sec%direction(jseg) .NE. sec%direction(1)) iptegrid = 1   
528        jseg = jseg+1
529     ENDDO 
530     IF(lk_mpp) CALL mpp_sum(iptegrid)
531
532     !---------------------------!
533     !  COMPUTE TRANSPORT        !
534     !---------------------------!
535     IF(sec%nb_point .NE. 0)THEN   
536
537        !----------------------------------------------------------------------------------------------------
538        !Compute sign for velocities:
539        !
540        !convention:
541        !   non horizontal section: sens + is toward left hand of section
542        !       horizontal section: sens + is toward north of section
543        !
544        !
545        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
546        !       ----------------      -----------------     ---------------             --------------
547        !
548        !   isgnv=1         sens +     
549        !  ______         _____             ______(i1,j1)      (i1,j1)                                     
550        !(i0,j0) |           //|            |                  |                   sens +   
551        !        | isgnu=1  // |            |isgnu=1           |isgnu=1              /|\
552        !        |_______  //         ______|    \\            | ---\                 |
553        !               |             | isgnv=-1  \\ |         | ---/ sens +       ____________(i1,j1)
554        !               |             |          __\\|         |               (i0,j0)
555        !               |(i1,j1)      |(i0,j0)   sens +        |                      isgnv=1                                 
556        !                                                      (i0,j0)
557        !----------------------------------------------------------------------------------------------------
558        isgnu = 1
559        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
560        ELSE                                ; isgnv =  1
561        ENDIF
562
563        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
564
565        !--------------------------------------!
566        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
567        !--------------------------------------!
568        DO jseg=1,MAX(sec%nb_point-1,0)
569             
570           !------------------------------------------
571           ! Select good coordinate to have velocity of the segment
572           !
573           !              |U CASE(3)
574           !      CASE(0) |                     j+1
575           !    ____V_____|F____V_CASE(1)__
576           !              |
577           !              |                     j
578           !      T       |U CASE(2)
579           !       i      |      i+1
580           !-------------------------------------------
581           SELECT CASE( sec%direction(jseg) )
582           CASE(0)  ;   k = sec%listPoint(jseg)
583           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
584           CASE(2)  ;   k = sec%listPoint(jseg)
585           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
586           END SELECT
587
588           !-------------------------------
589           !  LOOP ON THE DENSITY CLASSES |
590           !-------------------------------
591           !The computation is made for each density class
592           DO jclass=1,MAX(1,sec%nb_class-1)
593
594              ztransp1=0._wp ; zrhoi1=0._wp ; zrhop1=0._wp ; ztemp1=0._wp ;zsal1=0._wp
595              ztransp2=0._wp ; zrhoi2=0._wp ; zrhop2=0._wp ; ztemp2=0._wp ;zsal2=0._wp
596   
597              !---------------------------|
598              !     LOOP ON THE LEVEL     |
599              !---------------------------|
600              !Sum of the transport on the vertical
601              DO jk=1,jpk
602                   
603
604                 ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point
605                 SELECT CASE( sec%direction(jseg) )
606                 CASE(0,1)
607                    ztn   = interp(k%I,k%J,jk,'V',tn)
608                    zsn   = interp(k%I,k%J,jk,'V',sn)
609                    zrhop = interp(k%I,k%J,jk,'V',rhop)
610                    zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
611                    zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
612                 CASE(2,3)
613                    ztn   = interp(k%I,k%J,jk,'U',tn)
614                    zsn   = interp(k%I,k%J,jk,'U',sn)
615                    zrhop = interp(k%I,k%J,jk,'U',rhop)
616                    zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
617                    zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
618                 END SELECT
619
620                 zfsdep= gdept(k%I,k%J,jk)
621 
622                 !----------------------------------------------!
623                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
624                 !----------------------------------------------!
625           IF(ld_debug)write(narea+200,*)k ,jk
626 
627                 IF ( (    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.    &
628                           (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.    &
629                           ( sec%zsigp(jclass) .EQ. 99.)) .AND.                 &
630                           ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
631                           (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.    &
632                           ( sec%zsigi(jclass) .EQ. 99.)) .AND.                 &
633                           ((( zsn .GT. sec%zsal(jclass)) .AND.                &
634                           (   zsn .LE. sec%zsal(jclass+1))) .OR.              &
635                           ( sec%zsal(jclass) .EQ. 99.)) .AND.                 &
636                           ((( ztn .GE. sec%ztem(jclass)) .AND.                &
637                           (   ztn .LE. sec%ztem(jclass+1))) .OR.              &
638                           ( sec%ztem(jclass) .EQ.99.)) .AND.                  &
639                           ((( zfsdep .GE. sec%zlay(jclass)) .AND.            &
640                           (   zfsdep .LE. sec%zlay(jclass+1))) .OR.          &
641                           ( sec%zlay(jclass) .EQ. 99. ))))   THEN
642
643
644                    !compute velocity with the good sens
645                    SELECT CASE( sec%direction(jseg) )
646                    CASE(0,1) 
647                       zumid=0.
648                       zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
649                    CASE(2,3)
650                       zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
651                       zvmid=0.
652                    END SELECT
653
654                    !velocity* cell's lenght * cell's thickness
655                    zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &
656                           zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)
657
658#if ! defined key_vvl
659                    !add transport due to free surface
660                    IF( jk==1 )THEN
661                       zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
662                                         zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
663                    ENDIF
664#endif
665           !IF(ld_debug)write(narea+200,*)k ,jk,zTnorm             
666                    !COMPUTE TRANSPORT
667                    !zTnorm=transport through one cell for one class
668                    !ztransp1 or ztransp2=transport through one cell i
669                    !                     for one class for one direction
670                    IF( zTnorm .GE. 0 )THEN
671
672                       ztransp1=zTnorm+ztransp1
673 
674                       IF ( sec%llstrpond ) THEN
675                          ztemp1 = ztemp1  + zTnorm * ztn 
676                          zsal1  = zsal1   + zTnorm * zsn
677                          zrhoi1 = zrhoi1  + zTnorm * zrhoi
678                          zrhop1 = zrhop1  + zTnorm * zrhop
679                       ENDIF
680
681                    ELSE
682
683                       ztransp2=(zTnorm)+ztransp2
684
685                       IF ( sec%llstrpond ) THEN
686                          ztemp2 = ztemp2  + zTnorm * ztn 
687                          zsal2  = zsal2   + zTnorm * zsn
688                          zrhoi2 = zrhoi2  + zTnorm * zrhoi
689                          zrhop2 = zrhop2  + zTnorm * zrhop
690                       ENDIF
691                    ENDIF
692 
693           
694                 ENDIF ! end of density test
695              ENDDO!end of loop on the level
696
697              !ZSUM=TRANSPORT FOR EACH CLASSES FOR THE  DIRECTIONS
698              !---------------------------------------------------
699              zsum(1,jclass)     = zsum(1,jclass)+ztransp1
700              zsum(2,jclass)     = zsum(2,jclass)+ztransp2
701              IF( sec%llstrpond )THEN
702                 zsum(3 ,jclass) = zsum( 3,jclass)+zrhoi1
703                 zsum(4 ,jclass) = zsum( 4,jclass)+zrhoi2
704                 zsum(5 ,jclass) = zsum( 5,jclass)+zrhop1
705                 zsum(6 ,jclass) = zsum( 6,jclass)+zrhop2
706                 zsum(7 ,jclass) = zsum( 7,jclass)+ztemp1
707                 zsum(8 ,jclass) = zsum( 8,jclass)+ztemp2
708                 zsum(9 ,jclass) = zsum( 9,jclass)+zsal1
709                 zsum(10,jclass) = zsum(10,jclass)+zsal2
710              ENDIF
711   
712           ENDDO !end of loop on the density classes
713
714#if defined key_ice_lim
715
716           !ICE CASE   
717           !------------
718           IF( sec%ll_ice_section )THEN
719              SELECT CASE (sec%direction(jseg))
720              CASE(0)
721                 zumid_ice = 0
722                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
723              CASE(1)
724                 zumid_ice = 0
725                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
726              CASE(2)
727                 zvmid_ice = 0
728                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
729              CASE(3)
730                 zvmid_ice = 0
731                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
732              END SELECT
733   
734              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
735   
736              IF( zTnorm .GE. 0)THEN
737                 zice_vol_pos = (zTnorm)*   &
738                                      (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
739                                     *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
740                                       hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
741                                      +zice_vol_pos
742                 zice_surf_pos = (zTnorm)*   &
743                                       (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
744                                      +zice_surf_pos
745              ELSE
746                 zice_vol_neg=(zTnorm)*   &
747                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
748                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
749                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
750                                  +zice_vol_neg
751                 zice_surf_neg=(zTnorm)*   &
752                                    (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
753                                     +zice_surf_neg
754              ENDIF
755   
756              zsum(11,1) = zsum(11,1)+zice_vol_pos
757              zsum(12,1) = zsum(12,1)+zice_vol_neg
758              zsum(13,1) = zsum(13,1)+zice_surf_pos
759              zsum(14,1) = zsum(14,1)+zice_surf_neg
760   
761           ENDIF !end of ice case
762#endif
763 
764        ENDDO !end of loop on the segment
765
766
767     ELSE  !if nb_inmesh=0
768        zsum(1:2,:)=0.
769        IF (sec%llstrpond) zsum(3:10,:)=0.
770        zsum( 11:14,:)=0.
771     ENDIF   !end of nb_inmesh=0 case
772
773     !------------------!
774     !SUM ON ALL PROCS  !
775     !------------------!
776     IF( lk_mpp )THEN
777         ish(1)  =  nb_type_class*nb_class_max ;  ish2(1)=nb_type_class ; ish2(2)=nb_class_max
778         zwork(:)= RESHAPE(zsum(:,:), ish )
779         CALL mpp_sum(zwork, ish(1))
780         zsum(:,:)= RESHAPE(zwork,ish2)
781     ENDIF
782     
783     !-------------------------------|
784     !FINISH COMPUTING TRANSPORTS    |
785     !-------------------------------|
786     DO jclass=1,MAX(1,sec%nb_class-1)
787        sec%transport(1,jclass)=sec%transport(1,jclass)+zsum(1,jclass)*1.E-6
788        sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6
789        IF( sec%llstrpond ) THEN
790           IF( zsum(1,jclass) .NE. 0 ) THEN
791              sec%transport(3,jclass)=sec%transport(3,jclass)+zsum(3,jclass)/zsum(1,jclass)
792              sec%transport(5,jclass)=sec%transport(5,jclass)+zsum(5,jclass)/zsum(1,jclass)
793              sec%transport(7,jclass)=sec%transport(7,jclass)+zsum(7,jclass)
794              sec%transport(9,jclass)=sec%transport(9,jclass)+zsum(9,jclass)
795           ELSE
796              sec%transport(3,jclass)=0.
797              sec%transport(5,jclass)=0.
798              sec%transport(7,jclass)=0.
799              sec%transport(9,jclass)=0.
800           ENDIF
801           IF( zsum(2,jclass) .NE. 0 )THEN
802              sec%transport( 4,jclass)=sec%transport( 4,jclass)+zsum( 4,jclass)/zsum(2,jclass)
803              sec%transport( 6,jclass)=sec%transport( 6,jclass)+zsum( 6,jclass)/zsum(2,jclass)
804              sec%transport( 8,jclass)=sec%transport( 8,jclass)+zsum( 8,jclass)
805              sec%transport(10,jclass)=sec%transport(10,jclass)+zsum(10,jclass)
806           ELSE
807              sec%transport( 4,jclass)=0.
808              sec%transport( 6,jclass)=0.
809              sec%transport( 8,jclass)=0.
810              sec%transport(10,jclass)=0.
811           ENDIF
812        ELSE
813           sec%transport( 3,jclass)=0.
814           sec%transport( 4,jclass)=0.
815           sec%transport( 5,jclass)=0.
816           sec%transport( 6,jclass)=0.
817           sec%transport( 7,jclass)=0.
818           sec%transport( 8,jclass)=0.
819           sec%transport(10,jclass)=0.
820        ENDIF
821     ENDDO   
822
823     IF( sec%ll_ice_section ) THEN
824        sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6
825        sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6
826        sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6
827        sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6
828     ENDIF
829
830  END SUBROUTINE transport
831 
832  SUBROUTINE dia_dct_wri(kt,jsec,sec)
833     !!-------------------------------------------------------------
834     !! Write transport output in numdct
835     !!
836     !! Purpose: Write  transports in ascii files
837     !!
838     !! Method:
839     !!        1. Write volume transports in "volume_transport"
840     !!           Unit: Sv : Surface * Velocity / 1.e6
841     !!
842     !!        2. Write heat transports in "heat_transport"
843     !!           Unit: Peta W : Surface * Velocity * T * rhau * Cp / 1.e15
844     !!
845     !!        3. Write salt transports in "salt_transport"
846     !!           Unit: Mega m^3 / s : Surface * Velocity * S / 1.e6
847     !!
848     !!-------------------------------------------------------------
849     !!arguments
850     INTEGER, INTENT(IN)          :: kt
851     TYPE(SECTION), INTENT(INOUT) :: sec
852     INTEGER ,INTENT(IN)          :: jsec
853
854     !!local declarations
855     REAL(wp) ,DIMENSION(nb_type_class):: zsumclass
856     INTEGER               :: jcl,ji        !Dummy loop
857     CHARACTER(len=2)      :: classe
858     REAL(wp)              :: zbnd1,zbnd2
859     !!-------------------------------------------------------------
860     
861     zsumclass(:)=0._wp
862               
863     DO jcl=1,MAX(1,sec%nb_class-1)
864
865        ! Mean computation
866        sec%transport(:,jcl)=sec%transport(:,jcl)/(nn_dctwri/nn_dct)
867        classe   = 'N       '
868        zbnd1   = 0._wp
869        zbnd2   = 0._wp
870        zsumclass(1:nb_type_class)=zsumclass(1:nb_type_class)+sec%transport(1:nb_type_class,jcl)
871
872   
873        !insitu density classes transports
874        IF( ( sec%zsigi(jcl)   .NE. 99._wp ) .AND. &
875            ( sec%zsigi(jcl+1) .NE. 99._wp )       )THEN
876           classe = 'DI       '
877           zbnd1 = sec%zsigi(jcl)
878           zbnd2 = sec%zsigi(jcl+1)
879        ENDIF
880        !potential density classes transports
881        IF( ( sec%zsigp(jcl)   .NE. 99._wp ) .AND. &
882            ( sec%zsigp(jcl+1) .NE. 99._wp )       )THEN
883           classe = 'DP      '
884           zbnd1 = sec%zsigp(jcl)
885           zbnd2 = sec%zsigp(jcl+1)
886        ENDIF
887        !depth classes transports
888        IF( ( sec%zlay(jcl)    .NE. 99._wp ) .AND. &
889            ( sec%zlay(jcl+1)  .NE. 99._wp )       )THEN
890           classe = 'Z       '
891           zbnd1 = sec%zlay(jcl)
892           zbnd2 = sec%zlay(jcl+1)
893        ENDIF
894        !salinity classes transports
895        IF( ( sec%zsal(jcl) .NE. 99._wp    ) .AND. &
896            ( sec%zsal(jcl+1) .NE. 99._wp  )       )THEN
897           classe = 'S       '
898           zbnd1 = sec%zsal(jcl)
899           zbnd2 = sec%zsal(jcl+1)   
900        ENDIF
901        !temperature classes transports
902        IF( ( sec%ztem(jcl+1) .NE. 99.     ) .AND. &
903            ( sec%ztem(jcl+1) .NE. 99.     )       ) THEN
904           classe = 'T       '
905           zbnd1 = sec%ztem(jcl)
906           zbnd2 = sec%ztem(jcl+1)
907        ENDIF
908                 
909        !write volume transport per class
910        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, &
911                              jcl,classe,zbnd1,zbnd2,&
912                              sec%transport(1,jcl),sec%transport(2,jcl), &
913                              sec%transport(1,jcl)+sec%transport(2,jcl)
914
915        IF( sec%llstrpond )THEN
916
917           !write heat transport per class:
918           WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection,  &
919                              jcl,classe,zbnd1,zbnd2,&
920                              sec%transport(7,jcl)*1000._wp*rcp/1.e15,sec%transport(8,jcl)*1000._wp*rcp/1.e15, &
921                              ( sec%transport(7,jcl)+sec%transport(8,jcl) )*1000._wp*rcp/1.e15
922           !write salt transport per class
923           WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection,  &
924                              jcl,classe,zbnd1,zbnd2,&
925                              sec%transport(9,jcl)*1000._wp/1.e9,sec%transport(10,jcl)*1000._wp/1.e9,&
926                              (sec%transport(9,jcl)+sec%transport(10,jcl))*1000._wp/1.e9
927        ENDIF
928
929     ENDDO
930
931     zbnd1 = 0._wp
932     zbnd2 = 0._wp
933     jcl=0
934
935     !write total volume transport
936     WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, &
937                           jcl,"total",zbnd1,zbnd2,&
938                           zsumclass(1),zsumclass(2),zsumclass(1)+zsumclass(2)
939
940     IF( sec%llstrpond )THEN
941
942        !write total heat transport
943        WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
944                           jcl,"total",zbnd1,zbnd2,&
945                           zsumclass(7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,(zsumclass(7)+zsumclass(8) )* 1000._wp*rcp/1.e15
946        !write total salt transport
947        WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
948                           jcl,"total",zbnd1,zbnd2,&
949                           zsumclass(9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,(zsumclass(9)+zsumclass(10))*1000._wp/1.e9
950     ENDIF
951
952     
953     IF ( sec%ll_ice_section) THEN
954        !write total ice volume transport
955        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,&
956                              jcl,"ice_vol",zbnd1,zbnd2,&
957                              sec%transport(9,1),sec%transport(10,1),&
958                              sec%transport(9,1)+sec%transport(10,1)
959        !write total ice surface transport
960        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,&
961                              jcl,"ice_surf",zbnd1,zbnd2,&
962                              sec%transport(11,1),sec%transport(12,1), &
963                              sec%transport(11,1)+sec%transport(12,1) 
964     ENDIF
965                                             
966118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
967119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
968
969  END SUBROUTINE dia_dct_wri
970
971  FUNCTION interp(ki, kj, kk, cd_point, ptab)
972  !!----------------------------------------------------------------------
973  !!
974  !!   Purpose: compute Temperature/Salinity/density on U-point or V-point
975  !!
976  !!    |    I       |    I+1     |
977  !!    |            |            |
978  !!  --------------------------------     1.Compute Tbis:
979  !!    |            |            |          interpolate T(I,J,K) et T(I,J,K+1)
980  !!    |            |            |
981  !!    |            |            |        2. 1.Compute Temp/Salt/Rho à U point:
982  !! K  |    T       |            |          interpolate Tbis et T(I+U,J,K+1) 
983  !!    |            |            |
984  !!    |            |            |
985  !!    |            |            |
986  !!  --------------------------------
987  !!    |            |            |
988  !!    |            |            |
989  !!    |            |            |
990  !!K+1 |    Tbis    U     T      |
991  !!    |            |            |
992  !!    |    T       |            |
993  !!    |            |------------|
994  !!    |            | partials   |
995  !!    |            |  steps     |
996  !!  --------------------------------
997  !!
998  !!----------------------------------------------------------------------
999  !*arguments
1000  INTEGER, INTENT(IN)                          :: ki, kj, kk
1001  CHARACTER(len=1), INTENT(IN)                 :: cd_point
1002  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab
1003  REAL(wp)                                     :: interp
1004
1005  !*local declations
1006  INTEGER :: ii1, ij1, ii2, ij2
1007  REAL(wp):: ze3w, zfse3, zmax1, zmax2, zbis
1008  !!----------------------------------------------------------------------
1009
1010  IF( cd_point=='U' )THEN
1011     ii1 = ki    ; ij1 = kj 
1012     ii2 = ki+1  ; ij2 = kj 
1013  ELSE ! cd_point=='V'
1014     ii1 = ki    ; ij1 = kj 
1015     ii2 = ki    ; ij2 = kj+1 
1016  ENDIF
1017
1018#if defined key_vvl
1019
1020  ze3w  = fsve3w(ii2,ij2,kk) - fsve3w(ii1,ij1,kk)
1021
1022  zfse3 = fsve3w(ii1,ij1,kk) * ( 1 + sshn(ii2,ij2) * mut(ii2,ij2,kk) )
1023  zmax1 = ( fse3w(ii2,ij2,kk) - zfse3 )/ fse3w(ii2,ij2,kk)
1024 
1025  zfse3 = fsve3w(ii2,ij2,kk) * ( 1 + sshn(ii1,ij1) * mut(ii1,ij1,kk) )
1026  zmax2 = ( fse3w(ii1,ij1,kk) - zfse3 )/ fse3w(ii1,ij1,kk)
1027
1028#else
1029
1030  ze3w  = fse3w(ii2,ij2,kk)-fse3w(ii1,ij1,kk) 
1031  zmax1 =  ze3w / fse3w(ii2,ij2,kk)
1032  zmax2 = -ze3w / fse3w(ii1,ij1,kk)
1033
1034#endif
1035
1036  IF(kk .NE. 1)THEN
1037
1038     IF( ze3w >= 0. )THEN 
1039        !zbis
1040        zbis = ptab(ii2,ij2,kk) + zmax1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
1041        ! result
1042        interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + zbis )
1043     ELSE
1044        !zbis
1045        zbis = ptab(ii1,ij1,kk) + zmax2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1046        ! result
1047        interp = umask(ii1,ij1,kk) * 0.5*( zbis +  ptab(ii2,ij2,kk) )
1048     ENDIF   
1049
1050  ELSE
1051     interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + ptab(ii2,ij2,kk) )
1052  ENDIF
1053
1054  END FUNCTION interp
1055
1056#else
1057   !!----------------------------------------------------------------------
1058   !!   Default option :                                       Dummy module
1059   !!----------------------------------------------------------------------
1060   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
1061   PUBLIC
1062CONTAINS
1063
1064   SUBROUTINE dia_dct_init          ! Dummy routine
1065      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
1066   END SUBROUTINE dia_dct_init
1067
1068   SUBROUTINE dia_dct( kt )           ! Dummy routine
1069      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
1070      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
1071   END SUBROUTINE dia_dct
1072#endif
1073
1074END MODULE diadct
Note: See TracBrowser for help on using the repository browser.