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

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

add routines for online transports computing

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