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

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

modify formula of heat and salt transport

  • Property svn:executable set to *
File size: 46.5 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              ENDIF
338
339           ENDDO
340     
341           secs(jsec)%nb_point=iptloc !store number of section's points
342
343           !debug
344           !-----
345           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
346              WRITE(numout,*)"      List of points selected by the proc:"
347              DO jpt = 1,iptloc
348                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
349                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
350                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
351              ENDDO
352           ENDIF
353
354           !remove redundant points between processors
355           !look NEMO documentation for more explanations
356!           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
357           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1)  ) lldebug = .TRUE.
358           IF( iptloc .NE. 0 )THEN
359              CALL removepoints(secs(jsec),'I','top_list',lldebug)
360              CALL removepoints(secs(jsec),'I','bot_list',lldebug)
361              CALL removepoints(secs(jsec),'J','top_list',lldebug)
362              CALL removepoints(secs(jsec),'J','bot_list',lldebug)
363           ENDIF
364
365           !debug
366           !-----
367           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
368              WRITE(numout,*)"      List of points after removepoints:"
369              DO jpt = 1,iptloc
370                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
371                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
372                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
373              ENDDO
374           ENDIF
375
376        ELSE  ! iptglo = 0
377           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )&
378              WRITE(numout,*)'   Not points for this section.'
379        ENDIF
380
381     ENDDO !end of the loop on jsec
382 
383     nb_sec = jsec-1   !number of section read in the file
384
385  END SUBROUTINE readsec
386
387  SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
388     !!---------------------------------------------------------------------------
389     !!             *** function removepoints
390     !!
391     !!   ** Purpose ::
392     !!              remove points which are common to 2 procs
393     !!
394     !!
395     !----------------------------------------------------------------------------
396     !! * arguments
397     TYPE(SECTION),INTENT(INOUT) :: sec
398     CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
399     CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
400     LOGICAL,INTENT(IN)          :: ld_debug                     
401
402     !! * Local variables
403     INTEGER :: iextr         ,& !extremity of listpoint that we verify
404                iind          ,& !coord     of listpoint that we verify
405                itest         ,& !indice value of the side of the domain
406                                 !where points could be redundant
407                isgn          ,& !way of course in listpoint
408                istart,iend      !first and last points selected in listpoint
409     INTEGER :: jpoint   =0      !loop on list points
410     INTEGER,DIMENSION(nb_point_max)   :: idirec !contains temporare sec%direction
411     INTEGER,DIMENSION(2,nb_point_max) :: icoord !contains temporare sec%listpoint
412     !----------------------------------------------------------------------------
413     IF( ld_debug )WRITE(numout,*)'      -------------------------'
414     IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
415
416     !iextr=extremity of list_point that we verify
417     IF      ( cdextr=='bot_list' )THEN ; iextr=1            ; isgn=1
418     ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
419     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdextr")
420     ENDIF
421 
422     !which coordinate shall we verify ?
423     IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1
424     ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2
425     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind") 
426     ENDIF
427
428     IF( ld_debug )THEN
429        WRITE(numout,*)'      case: coord/list extr/domain side'
430        WRITE(numout,*)'      ', cdind,' ',cdextr,' ',itest
431        WRITE(numout,*)'      Actual number of points: ',sec%nb_point
432     ENDIF
433
434     icoord(1,1:nb_point_max) = sec%listPoint%I
435     icoord(2,1:nb_point_max) = sec%listPoint%J
436     idirec                   = sec%direction
437     sec%listPoint            = POINT_SECTION(0,0)
438     sec%direction            = 0
439
440
441     jpoint=iextr+isgn
442     DO WHILE ( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point .AND. &
443             icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest)
444             jpoint=jpoint+isgn
445     ENDDO
446     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
447     ELSE                        ; istart=1        ; iend=jpoint+1
448     ENDIF
449     sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
450     sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
451     sec%direction(1:1+iend-istart)   = idirec(istart:iend)
452     sec%nb_point                     = iend-istart+1
453     
454     IF( ld_debug )THEN
455        WRITE(numout,*)'      Number of points after removepoints :',sec%nb_point
456        WRITE(numout,*)'      sec%direction after removepoints :',sec%direction(1:sec%nb_point)
457     ENDIF
458
459  END SUBROUTINE removepoints
460
461  SUBROUTINE transport(sec,ld_debug)
462     !!---------------------------------------------------------------------
463     !!                     ***  ROUTINE transport  ***
464     !!
465     !!  ** Purpose : Compute the transport trough a sec
466     !!
467     !!  ** Method  :Transport through a given section is equal to the sum of transports
468     !!              computed on each proc.
469     !!              On each proc,transport is equal to the sum of transport computed through
470     !!               segments linking each points of sec%listPoint  with the next one.           
471     !!              There are several loops:                 
472     !!              loop on the density/temperature/salinity/level classes
473     !!              loop on the segment between 2 nodes
474     !!              loop on the level jk
475     !!              test on the density/temperature/salinity/level
476     !!
477     !! ** Output: sec%transport: transport in the 2 direction and temperature,
478     !!                       salinity, density meaned by the transport
479     !!
480     !!
481     !!-------------------------------------------------------------------
482     !! * Arguments
483     TYPE(SECTION),INTENT(INOUT) :: sec
484     LOGICAL      ,INTENT(IN)    :: ld_debug
485   
486     !! * Local variables
487     INTEGER             :: jk,jseg,jclass,   &!loop on level/segment/classes
488                            iptegrid,         &!=0 if section is constant along i/j
489                            isgnu  , isgnv     !
490     INTEGER :: ii, ij ! local integer
491     REAL(wp):: zumid        , zvmid        ,&!U/V velocity on a cell segment
492                zumid_ice    , zvmid_ice    ,&!U/V ice velocity
493                zTnorm                      ,&!transport of velocity through one cell's sides
494                ztransp1     , ztransp2     ,&!total        transport in sens 1 and 2
495                ztemp1       , ztemp2       ,&!temperature  transport     "
496                zrhoi1       , zrhoi2       ,&!mass         transport     "
497                zrhop1       , zrhop2       ,&!mass         transport     "
498                zsal1        , zsal2        ,&!salinity     transport     "
499                zice_vol_pos , zice_vol_neg ,&!volumic  ice transport     "
500                zice_surf_pos, zice_surf_neg  !surfacic ice transport     "
501     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
502     REAL(wp) :: aj,i0,j0,i1,j1,i,j
503
504     TYPE(POINT_SECTION) :: k
505     INTEGER,DIMENSION(1):: ish
506     INTEGER,DIMENSION(2):: ish2
507     REAL(wp),DIMENSION(nb_type_class,nb_class_max)::zsum
508     REAL(wp),DIMENSION(nb_type_class*nb_class_max)::zwork  ! temporary vector for mpp_sum
509     !!--------------------------------------------------------
510
511     IF( ld_debug )WRITE(numout,*)'      Compute transport'
512
513     !----------------!
514     ! INITIALIZATION !
515     !----------------!
516     zsum    = 0._wp
517     zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp
518     zice_vol_pos  = 0._wp ; zice_vol_neg  = 0._wp
519
520     !compute iptegrid: iptegrid=0 if section is constant along i or j
521     iptegrid = 0
522     jseg     = 1
523     DO WHILE( (jseg .LT. MAX(sec%nb_point-1,0)) .AND. iptegrid .NE. 1 )
524        IF(sec%direction(jseg) .NE. sec%direction(1)) iptegrid = 1   
525        jseg = jseg+1
526     ENDDO 
527     IF(lk_mpp) CALL mpp_sum(iptegrid)
528
529     !---------------------------!
530     !  COMPUTE TRANSPORT        !
531     !---------------------------!
532     IF(sec%nb_point .NE. 0)THEN   
533
534        !----------------------------------------------------------------------------------------------------
535        !Compute sign for velocities:
536        !
537        !convention:
538        !   non horizontal section: sens + is toward left hand of section
539        !       horizontal section: sens + is toward north of section
540        !
541        !
542        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
543        !       ----------------      -----------------     ---------------             --------------
544        !
545        !   isgnv=1         sens +     
546        !  ______         _____             ______(i1,j1)      (i1,j1)                                     
547        !(i0,j0) |           //|            |                  |                   sens +   
548        !        | isgnu=1  // |            |isgnu=1           |isgnu=1              /|\
549        !        |_______  //         ______|    \\            | ---\                 |
550        !               |             | isgnv=-1  \\ |         | ---/ sens +       ____________(i1,j1)
551        !               |             |          __\\|         |               (i0,j0)
552        !               |(i1,j1)      |(i0,j0)   sens +        |                      isgnv=1                                 
553        !                                                      (i0,j0)
554        !----------------------------------------------------------------------------------------------------
555        isgnu = 1
556        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
557        ELSE                                ; isgnv =  1
558        ENDIF
559
560        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
561
562        !--------------------------------------!
563        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
564        !--------------------------------------!
565        DO jseg=1,MAX(sec%nb_point-1,0)
566             
567           !------------------------------------------
568           ! Select good coordinate to have velocity of the segment
569           !
570           !              |U CASE(3)
571           !      CASE(0) |                     j+1
572           !    ____V_____|F____V_CASE(1)__
573           !              |
574           !              |                     j
575           !      T       |U CASE(2)
576           !       i      |      i+1
577           !-------------------------------------------
578           SELECT CASE( sec%direction(jseg) )
579           CASE(0)  ;   k = sec%listPoint(jseg)
580           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
581           CASE(2)  ;   k = sec%listPoint(jseg)
582           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
583           END SELECT
584             
585           !-------------------------------
586           !  LOOP ON THE DENSITY CLASSES |
587           !-------------------------------
588           !The computation is made for each density class
589           DO jclass=1,MAX(1,sec%nb_class-1)
590
591              ztransp1=0._wp ; zrhoi1=0._wp ; zrhop1=0._wp ; ztemp1=0._wp ;zsal1=0._wp
592              ztransp2=0._wp ; zrhoi2=0._wp ; zrhop2=0._wp ; ztemp2=0._wp ;zsal2=0._wp
593   
594              !---------------------------|
595              !     LOOP ON THE LEVEL     |
596              !---------------------------|
597              !Sum of the transport on the vertical
598              DO jk=1,jpk
599                   
600
601                 ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point
602                 SELECT CASE( sec%direction(jseg) )
603                 CASE(0,1)
604                    ztn   = interp(k%I,k%J,jk,'V',tn)
605                    zsn   = interp(k%I,k%J,jk,'V',sn)
606                    zrhop = interp(k%I,k%J,jk,'V',rhop)
607                    zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
608                    zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
609                 CASE(2,3)
610                    ztn   = interp(k%I,k%J,jk,'U',tn)
611                    zsn   = interp(k%I,k%J,jk,'U',sn)
612                    zrhop = interp(k%I,k%J,jk,'U',rhop)
613                    zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
614                    zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
615                 END SELECT
616
617                 zfsdep= gdept(k%I,k%J,jk)
618 
619                 !----------------------------------------------!
620                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
621                 !----------------------------------------------!
622 
623                 IF ( (    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.    &
624                           (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.    &
625                           ( sec%zsigp(jclass) .EQ. 99.)) .AND.                 &
626                           ((( zrhoi .GE. (sec%zsigi(jclass)  )) .AND.    &
627                           (   zrhoi .LE. (sec%zsigi(jclass+1)))) .OR.    &
628                           ( sec%zsigi(jclass) .EQ. 99.)) .AND.                 &
629                           ((( zsn .GT. sec%zsal(jclass)) .AND.                &
630                           (   zsn .LE. sec%zsal(jclass+1))) .OR.              &
631                           ( sec%zsal(jclass) .EQ. 99.)) .AND.                 &
632                           ((( ztn .GE. sec%ztem(jclass)) .AND.                &
633                           (   ztn .LE. sec%ztem(jclass+1))) .OR.              &
634                           ( sec%ztem(jclass) .EQ.99.)) .AND.                  &
635                           ((( zfsdep .GE. sec%zlay(jclass)) .AND.            &
636                           (   zfsdep .LE. sec%zlay(jclass+1))) .OR.          &
637                           ( sec%zlay(jclass) .EQ. 99. ))))   THEN
638
639
640                    !compute velocity with the good sens
641                    SELECT CASE( sec%direction(jseg) )
642                    CASE(0,1) 
643                       zumid=0.
644                       zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
645                    CASE(2,3)
646                       zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
647                       zvmid=0.
648                    END SELECT
649
650                    !velocity* cell's lenght * cell's thickness
651                    zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &
652                           zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)
653
654#if ! defined key_vvl
655                    !add transport due to free surface
656                    IF( jk==1 )THEN
657                       zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
658                                         zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
659                    ENDIF
660#endif
661                    !COMPUTE TRANSPORT
662                    !zTnorm=transport through one cell for one class
663                    !ztransp1 or ztransp2=transport through one cell i
664                    !                     for one class for one direction
665                    IF( zTnorm .GE. 0 )THEN
666
667                       ztransp1=zTnorm+ztransp1
668 
669                       IF ( sec%llstrpond ) THEN
670                          ztemp1 = ztemp1  + zTnorm * ztn 
671                          zsal1  = zsal1   + zTnorm * zsn
672                          zrhoi1 = zrhoi1  + zTnorm * zrhoi
673                          zrhop1 = zrhop1  + zTnorm * zrhop
674                       ENDIF
675
676                    ELSE
677
678                       ztransp2=(zTnorm)+ztransp2
679
680                       IF ( sec%llstrpond ) THEN
681                          ztemp2 = ztemp2  + zTnorm * ztn 
682                          zsal2  = zsal2   + zTnorm * zsn
683                          zrhoi2 = zrhoi2  + zTnorm * zrhoi
684                          zrhop2 = zrhop2  + zTnorm * zrhop
685                       ENDIF
686                    ENDIF
687 
688           
689                 ENDIF ! end of density test
690              ENDDO!end of loop on the level
691
692              !ZSUM=TRANSPORT FOR EACH CLASSES FOR THE  DIRECTIONS
693              !---------------------------------------------------
694              zsum(1,jclass)     = zsum(1,jclass)+ztransp1
695              zsum(2,jclass)     = zsum(2,jclass)+ztransp2
696              IF( sec%llstrpond )THEN
697                 zsum(3 ,jclass) = zsum( 3,jclass)+zrhoi1
698                 zsum(4 ,jclass) = zsum( 4,jclass)+zrhoi2
699                 zsum(5 ,jclass) = zsum( 5,jclass)+zrhop1
700                 zsum(6 ,jclass) = zsum( 6,jclass)+zrhop2
701                 zsum(7 ,jclass) = zsum( 7,jclass)+ztemp1
702                 zsum(8 ,jclass) = zsum( 8,jclass)+ztemp2
703                 zsum(9 ,jclass) = zsum( 9,jclass)+zsal1
704                 zsum(10,jclass) = zsum(10,jclass)+zsal2
705              ENDIF
706   
707           ENDDO !end of loop on the density classes
708
709#if defined key_ice_lim
710
711           !ICE CASE   
712           !------------
713           IF( sec%ll_ice_section )THEN
714              SELECT CASE (sec%direction(jseg))
715              CASE(0)
716                 zumid_ice = 0
717                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
718              CASE(1)
719                 zumid_ice = 0
720                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
721              CASE(2)
722                 zvmid_ice = 0
723                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
724              CASE(3)
725                 zvmid_ice = 0
726                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
727              END SELECT
728   
729              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
730   
731              IF( zTnorm .GE. 0)THEN
732                 zice_vol_pos = (zTnorm)*   &
733                                      (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
734                                     *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
735                                       hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
736                                      +zice_vol_pos
737                 zice_surf_pos = (zTnorm)*   &
738                                       (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
739                                      +zice_surf_pos
740              ELSE
741                 zice_vol_neg=(zTnorm)*   &
742                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
743                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
744                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
745                                  +zice_vol_neg
746                 zice_surf_neg=(zTnorm)*   &
747                                    (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
748                                     +zice_surf_neg
749              ENDIF
750   
751              zsum(11,1) = zsum(11,1)+zice_vol_pos
752              zsum(12,1) = zsum(12,1)+zice_vol_neg
753              zsum(13,1) = zsum(13,1)+zice_surf_pos
754              zsum(14,1) = zsum(14,1)+zice_surf_neg
755   
756           ENDIF !end of ice case
757#endif
758 
759        ENDDO !end of loop on the segment
760
761
762     ELSE  !if nb_inmesh=0
763        zsum(1:2,:)=0.
764        IF (sec%llstrpond) zsum(3:10,:)=0.
765        zsum( 11:14,:)=0.
766     ENDIF   !end of nb_inmesh=0 case
767
768     !------------------!
769     !SUM ON ALL PROCS  !
770     !------------------!
771     IF( lk_mpp )THEN
772         ish(1)  =  nb_type_class*nb_class_max ;  ish2(1)=nb_type_class ; ish2(2)=nb_class_max
773         zwork(:)= RESHAPE(zsum(:,:), ish )
774         CALL mpp_sum(zwork, ish(1))
775         zsum(:,:)= RESHAPE(zwork,ish2)
776     ENDIF
777     
778     !-------------------------------|
779     !FINISH COMPUTING TRANSPORTS    |
780     !-------------------------------|
781     DO jclass=1,MAX(1,sec%nb_class-1)
782        sec%transport(1,jclass)=sec%transport(1,jclass)+zsum(1,jclass)*1.E-6
783        sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6
784        IF( sec%llstrpond ) THEN
785           IF( zsum(1,jclass) .NE. 0 ) THEN
786              sec%transport(3,jclass)=sec%transport(3,jclass)+zsum(3,jclass)/zsum(1,jclass)
787              sec%transport(5,jclass)=sec%transport(5,jclass)+zsum(5,jclass)/zsum(1,jclass)
788              sec%transport(7,jclass)=sec%transport(7,jclass)+zsum(7,jclass)
789              sec%transport(9,jclass)=sec%transport(9,jclass)+zsum(9,jclass)
790           ELSE
791              sec%transport(3,jclass)=0.
792              sec%transport(5,jclass)=0.
793              sec%transport(7,jclass)=0.
794              sec%transport(9,jclass)=0.
795           ENDIF
796           IF( zsum(2,jclass) .NE. 0 )THEN
797              sec%transport( 4,jclass)=sec%transport( 4,jclass)+zsum( 4,jclass)/zsum(2,jclass)
798              sec%transport( 6,jclass)=sec%transport( 6,jclass)+zsum( 6,jclass)/zsum(2,jclass)
799              sec%transport( 8,jclass)=sec%transport( 8,jclass)+zsum( 8,jclass)
800              sec%transport(10,jclass)=sec%transport(10,jclass)+zsum(10,jclass)
801           ELSE
802              sec%transport( 4,jclass)=0.
803              sec%transport( 6,jclass)=0.
804              sec%transport( 8,jclass)=0.
805              sec%transport(10,jclass)=0.
806           ENDIF
807        ELSE
808           sec%transport( 3,jclass)=0.
809           sec%transport( 4,jclass)=0.
810           sec%transport( 5,jclass)=0.
811           sec%transport( 6,jclass)=0.
812           sec%transport( 7,jclass)=0.
813           sec%transport( 8,jclass)=0.
814           sec%transport(10,jclass)=0.
815        ENDIF
816     ENDDO   
817
818     IF( sec%ll_ice_section ) THEN
819        sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6
820        sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6
821        sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6
822        sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6
823     ENDIF
824
825  END SUBROUTINE transport
826 
827  SUBROUTINE dia_dct_wri(kt,jsec,sec)
828     !!-------------------------------------------------------------
829     !! Write transport output in numdct
830     !!
831     !! Purpose: Write  transports in ascii files
832     !!
833     !! Method:
834     !!        1. Write volume transports in "volume_transport"
835     !!           Unit: Sv : Surface * Velocity / 1.e6
836     !!
837     !!        2. Write heat transports in "heat_transport"
838     !!           Unit: Peta W : Surface * Velocity * T * rhau * Cp / 1.e15
839     !!
840     !!        3. Write salt transports in "salt_transport"
841     !!           Unit: Mega m^3 / s : Surface * Velocity * S / 1.e6
842     !!
843     !!-------------------------------------------------------------
844     !!arguments
845     INTEGER, INTENT(IN)          :: kt
846     TYPE(SECTION), INTENT(INOUT) :: sec
847     INTEGER ,INTENT(IN)          :: jsec
848
849     !!local declarations
850     REAL(wp) ,DIMENSION(nb_type_class):: zsumclass
851     INTEGER               :: jcl,ji        !Dummy loop
852     CHARACTER(len=2)      :: classe
853     REAL(wp)              :: zbnd1,zbnd2
854     !!-------------------------------------------------------------
855     
856     zsumclass(:)=0._wp
857               
858     DO jcl=1,MAX(1,sec%nb_class-1)
859
860        ! Mean computation
861        sec%transport(:,jcl)=sec%transport(:,jcl)/(nn_dctwri/nn_dct)
862        classe   = 'N       '
863        zbnd1   = 0._wp
864        zbnd2   = 0._wp
865        zsumclass(1:nb_type_class)=zsumclass(1:nb_type_class)+sec%transport(1:nb_type_class,jcl)
866
867   
868        !insitu density classes transports
869        IF( ( sec%zsigi(jcl)   .NE. 99._wp ) .AND. &
870            ( sec%zsigi(jcl+1) .NE. 99._wp )       )THEN
871           classe = 'DI       '
872           zbnd1 = sec%zsigi(jcl)
873           zbnd2 = sec%zsigi(jcl+1)
874        ENDIF
875        !potential density classes transports
876        IF( ( sec%zsigp(jcl)   .NE. 99._wp ) .AND. &
877            ( sec%zsigp(jcl+1) .NE. 99._wp )       )THEN
878           classe = 'DP      '
879           zbnd1 = sec%zsigp(jcl)
880           zbnd2 = sec%zsigp(jcl+1)
881        ENDIF
882        !depth classes transports
883        IF( ( sec%zlay(jcl)    .NE. 99._wp ) .AND. &
884            ( sec%zlay(jcl+1)  .NE. 99._wp )       )THEN
885           classe = 'Z       '
886           zbnd1 = sec%zlay(jcl)
887           zbnd2 = sec%zlay(jcl+1)
888        ENDIF
889        !salinity classes transports
890        IF( ( sec%zsal(jcl) .NE. 99._wp    ) .AND. &
891            ( sec%zsal(jcl+1) .NE. 99._wp  )       )THEN
892           classe = 'S       '
893           zbnd1 = sec%zsal(jcl)
894           zbnd2 = sec%zsal(jcl+1)   
895        ENDIF
896        !temperature classes transports
897        IF( ( sec%ztem(jcl+1) .NE. 99.     ) .AND. &
898            ( sec%ztem(jcl+1) .NE. 99.     )       ) THEN
899           classe = 'T       '
900           zbnd1 = sec%ztem(jcl)
901           zbnd2 = sec%ztem(jcl+1)
902        ENDIF
903                 
904        !write volume transport per class
905        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, &
906                              jcl,classe,zbnd1,zbnd2,&
907                              sec%transport(1,jcl),sec%transport(2,jcl), &
908                              sec%transport(1,jcl)+sec%transport(2,jcl)
909
910        IF( sec%llstrpond )THEN
911
912           !write heat transport per class:
913           WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection,  &
914                              jcl,classe,zbnd1,zbnd2,&
915                              sec%transport(7,jcl)*1000._wp*rcp/1.e15,sec%transport(8,jcl)*1000._wp*rcp/1.e15, &
916                              ( sec%transport(7,jcl)+sec%transport(8,jcl) )*1000._wp*rcp/1.e15
917           !write salt transport per class
918           WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection,  &
919                              jcl,classe,zbnd1,zbnd2,&
920                              sec%transport(9,jcl)*1000._wp/1.e9,sec%transport(10,jcl)*1000._wp/1.e9,&
921                              (sec%transport(9,jcl)+sec%transport(10,jcl))*1000._wp/1.e9
922        ENDIF
923
924     ENDDO
925
926     zbnd1 = 0._wp
927     zbnd2 = 0._wp
928     jcl=0
929
930     !write total volume transport
931     WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, &
932                           jcl,"total",zbnd1,zbnd2,&
933                           zsumclass(1),zsumclass(2),zsumclass(1)+zsumclass(2)
934
935     IF( sec%llstrpond )THEN
936
937        !write total heat transport
938        WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
939                           jcl,"total",zbnd1,zbnd2,&
940                           zsumclass(7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,(zsumclass(7)+zsumclass(8) )* 1000._wp*rcp/1.e15
941        !write total salt transport
942        WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
943                           jcl,"total",zbnd1,zbnd2,&
944                           zsumclass(9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,(zsumclass(9)+zsumclass(10))*1000._wp/1.e9
945     ENDIF
946
947     
948     IF ( sec%ll_ice_section) THEN
949        !write total ice volume transport
950        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,&
951                              jcl,"ice_vol",zbnd1,zbnd2,&
952                              sec%transport(9,1),sec%transport(10,1),&
953                              sec%transport(9,1)+sec%transport(10,1)
954        !write total ice surface transport
955        WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,&
956                              jcl,"ice_surf",zbnd1,zbnd2,&
957                              sec%transport(11,1),sec%transport(12,1), &
958                              sec%transport(11,1)+sec%transport(12,1) 
959     ENDIF
960                                             
961118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
962119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
963
964  END SUBROUTINE dia_dct_wri
965
966  FUNCTION interp(ki, kj, kk, cd_point, ptab)
967  !!----------------------------------------------------------------------
968  !!
969  !!   Purpose: compute Temperature/Salinity/density on U-point or V-point
970  !!
971  !!    |    I       |    I+1     |
972  !!    |            |            |
973  !!  --------------------------------     1.Compute Tbis:
974  !!    |            |            |          interpolate T(I,J,K) et T(I,J,K+1)
975  !!    |            |            |
976  !!    |            |            |        2. 1.Compute Temp/Salt/Rho à U point:
977  !! K  |    T       |            |          interpolate Tbis et T(I+U,J,K+1) 
978  !!    |            |            |
979  !!    |            |            |
980  !!    |            |            |
981  !!  --------------------------------
982  !!    |            |            |
983  !!    |            |            |
984  !!    |            |            |
985  !!K+1 |    Tbis    U     T      |
986  !!    |            |            |
987  !!    |    T       |            |
988  !!    |            |------------|
989  !!    |            | partials   |
990  !!    |            |  steps     |
991  !!  --------------------------------
992  !!
993  !!----------------------------------------------------------------------
994  !*arguments
995  INTEGER, INTENT(IN)                          :: ki, kj, kk
996  CHARACTER(len=1), INTENT(IN)                 :: cd_point
997  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab
998  REAL(wp)                                     :: interp
999
1000  !*local declations
1001  INTEGER :: ii1, ij1, ii2, ij2
1002  REAL(wp):: ze3w, zfse3, zmax1, zmax2, zbis
1003  !!----------------------------------------------------------------------
1004
1005  IF( cd_point=='U' )THEN
1006     ii1 = ki    ; ij1 = kj 
1007     ii2 = ki+1  ; ij2 = kj 
1008  ELSE ! cd_point=='V'
1009     ii1 = ki    ; ij1 = kj 
1010     ii2 = ki    ; ij2 = kj+1 
1011  ENDIF
1012
1013#if defined key_vvl
1014
1015  ze3w  = fsve3w(ii2,ij2,kk) - fsve3w(ii1,ij1,kk)
1016
1017  zfse3 = fsve3w(ii1,ij1,kk) * ( 1 + sshn(ii2,ij2) * mut(ii2,ij2,kk) )
1018  zmax1 = ( fse3w(ii2,ij2,kk) - zfse3 )/ fse3w(ii2,ij2,kk)
1019 
1020  zfse3 = fsve3w(ii2,ij2,kk) * ( 1 + sshn(ii1,ij1) * mut(ii1,ij1,kk) )
1021  zmax2 = ( fse3w(ii1,ij1,kk) - zfse3 )/ fse3w(ii1,ij1,kk)
1022
1023#else
1024
1025  ze3w  = fse3w(ii2,ij2,kk)-fse3w(ii1,ij1,kk) 
1026  zmax1 =  ze3w / fse3w(ii2,ij2,kk)
1027  zmax2 = -ze3w / fse3w(ii1,ij1,kk)
1028
1029#endif
1030
1031  IF(kk .NE. 1)THEN
1032
1033     IF( ze3w >= 0. )THEN 
1034        !zbis
1035        zbis = ptab(ii2,ij2,kk) + zmax1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
1036        ! result
1037        interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + zbis )
1038     ELSE
1039        !zbis
1040        zbis = ptab(ii1,ij1,kk) + zmax2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1041        ! result
1042        interp = umask(ii1,ij1,kk) * 0.5*( zbis +  ptab(ii2,ij2,kk) )
1043     ENDIF   
1044
1045  ELSE
1046     interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + ptab(ii2,ij2,kk) )
1047  ENDIF
1048
1049  END FUNCTION interp
1050
1051#else
1052   !!----------------------------------------------------------------------
1053   !!   Default option :                                       Dummy module
1054   !!----------------------------------------------------------------------
1055   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
1056CONTAINS
1057   SUBROUTINE dia_dct( kt )           ! Dummy routine
1058      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
1059      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
1060   END SUBROUTINE dia_dct
1061#endif
1062
1063END MODULE diadct
Note: See TracBrowser for help on using the repository browser.