source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 4 years ago

Merge of dev_CNRS_2017 into branch

  • Property svn:keywords set to Id
File size: 59.2 KB
Line 
1MODULE diadct
2   !!======================================================================
3   !!                       ***  MODULE  diadct  ***
4   !! Ocean diagnostics: Compute the transport trough a sec.
5   !!======================================================================
6   !! History :  OPA  ! 02/1999 (Y Drillet)  original code
7   !!                 ! 10/2001 (Y Drillet, R Bourdalle Badie)
8   !!   NEMO     1.0  ! 10/2005 (M Laborie) F90
9   !!            3.0  ! 04/2007 (G Garric) Ice sections
10   !!             -   ! 04/2007 (C Bricaud) test on sec%nb_point, initialisation of ztransp1,ztransp2,...
11   !!            3.4  ! 09/2011 (C Bricaud)
12   !!----------------------------------------------------------------------
13#if defined key_diadct
14   !!----------------------------------------------------------------------
15   !!   'key_diadct' :
16   !!----------------------------------------------------------------------
17   !!----------------------------------------------------------------------
18   !!   dia_dct      :  Compute the transport through a sec.
19   !!   dia_dct_init :  Read namelist.
20   !!   readsec      :  Read sections description and pathway
21   !!   removepoints :  Remove points which are common to 2 procs
22   !!   transport    :  Compute transport for each sections
23   !!   dia_dct_wri  :  Write tranports results in ascii files
24   !!   interp       :  Compute temperature/salinity/density at U-point or V-point
25   !!   
26   !!----------------------------------------------------------------------
27   USE oce             ! ocean dynamics and tracers
28   USE dom_oce         ! ocean space and time domain
29   USE phycst          ! physical constants
30   USE in_out_manager  ! I/O manager
31   USE daymod          ! calendar
32   USE dianam          ! build name of file
33   USE lib_mpp         ! distributed memory computing library
34#if defined key_lim3
35   USE ice
36#endif
37   USE domvvl
38   USE timing          ! preformance summary
39   USE wrk_nemo        ! working arrays
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC   dia_dct      ! routine called by step.F90
45   PUBLIC   dia_dct_init ! routine called by opa.F90
46   PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90
47   PRIVATE  readsec
48   PRIVATE  removepoints
49   PRIVATE  transport
50   PRIVATE  dia_dct_wri
51
52   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
53
54   INTEGER :: nn_dct        ! Frequency of computation
55   INTEGER :: nn_dctwri     ! Frequency of output
56   INTEGER :: nn_secdebug   ! 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 = 10
62   INTEGER, PARAMETER :: nb_3d_vars    = 3 
63   INTEGER, PARAMETER :: nb_2d_vars    = 2 
64   INTEGER            :: nb_sec 
65
66   TYPE POINT_SECTION
67      INTEGER :: I,J
68   END TYPE POINT_SECTION
69
70   TYPE COORD_SECTION
71      REAL(wp) :: lon,lat
72   END TYPE COORD_SECTION
73
74   TYPE SECTION
75      CHARACTER(len=60)                            :: name              ! name of the sec
76      LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
77                                                                       ! heat transports
78      LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
79      LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
80      TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
81      INTEGER                                      :: nb_class          ! number of boundaries for density classes
82      INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
83      CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
84      REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
85                                                      zsigp           ,&! potential density classes    (99 if you don't want)
86                                                      zsal            ,&! salinity classes   (99 if you don't want)
87                                                      ztem            ,&! temperature classes(99 if you don't want)
88                                                      zlay              ! level classes      (99 if you don't want)
89      REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output
90      REAL(wp)                                         :: slopeSection  ! slope of the section
91      INTEGER                                          :: nb_point      ! number of points in the section
92      TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
93   END TYPE SECTION
94
95   TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
96 
97   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d 
98   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d 
99
100   !!----------------------------------------------------------------------
101   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
102   !! $Id$
103   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
104   !!----------------------------------------------------------------------
105CONTAINS
106 
107  INTEGER FUNCTION diadct_alloc() 
108     !!----------------------------------------------------------------------
109     !!                   ***  FUNCTION diadct_alloc  ***
110     !!----------------------------------------------------------------------
111     INTEGER :: ierr(2) 
112     !!----------------------------------------------------------------------
113
114     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) ) 
115     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) ) 
116
117     diadct_alloc = MAXVAL( ierr ) 
118     IF( diadct_alloc /= 0 )   CALL ctl_warn('diadct_alloc: failed to allocate arrays') 
119 
120  END FUNCTION diadct_alloc 
121
122
123  SUBROUTINE dia_dct_init
124     !!---------------------------------------------------------------------
125     !!               ***  ROUTINE diadct  *** 
126     !!
127     !!  ** Purpose: Read the namelist parameters
128     !!              Open output files
129     !!
130     !!---------------------------------------------------------------------
131     NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
132     INTEGER  ::   ios                 ! Local integer output status for namelist read
133
134     IF( nn_timing == 1 )   CALL timing_start('dia_dct_init')
135
136     REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections
137     READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
138901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp )
139
140     REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
141     READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
142902  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp )
143     IF(lwm) WRITE ( numond, namdct )
144
145     IF( lwp ) THEN
146        WRITE(numout,*) " "
147        WRITE(numout,*) "diadct_init: compute transports through sections "
148        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
149        WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct
150        WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri
151
152        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
153                                            WRITE(numout,*)"       Debug section number: ", nn_secdebug 
154        ELSE IF ( nn_secdebug ==  0 )THEN ; WRITE(numout,*)"       No section to debug"
155        ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)"       Debug all sections"
156        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
157        ENDIF
158
159        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
160          &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
161
162     ENDIF
163
164     !Read section_ijglobal.diadct
165     CALL readsec
166
167     !open output file
168     IF( lwm ) THEN
169        CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
170        CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
171        CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
172     ENDIF
173
174     ! Initialise arrays to zero
175     transports_3d(:,:,:,:)=0.0 
176     transports_2d(:,:,:)  =0.0 
177
178     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init')
179     !
180  END SUBROUTINE dia_dct_init
181 
182 
183  SUBROUTINE dia_dct(kt)
184     !!---------------------------------------------------------------------
185     !!               ***  ROUTINE diadct  *** 
186     !!
187     !!  Purpose :: Compute section transports and write it in numdct files
188     !!   
189     !!  Method  :: All arrays initialised to zero in dct_init
190     !!             Each nn_dct time step call subroutine 'transports' for
191     !!               each section to sum the transports over each grid cell.
192     !!             Each nn_dctwri time step:
193     !!               Divide the arrays by the number of summations to gain
194     !!               an average value
195     !!               Call dia_dct_sum to sum relevant grid boxes to obtain
196     !!               totals for each class (density, depth, temp or sal)
197     !!               Call dia_dct_wri to write the transports into file
198     !!               Reinitialise all relevant arrays to zero
199     !!---------------------------------------------------------------------
200     INTEGER,INTENT(in)        ::kt
201     !
202     INTEGER             :: jsec,            &! loop on sections
203                            itotal            ! nb_sec_max*nb_type_class*nb_class_max
204     LOGICAL             :: lldebug =.FALSE.  ! debug a section 
205     
206     INTEGER , DIMENSION(1)             :: ish   ! tmp array for mpp_sum
207     INTEGER , DIMENSION(3)             :: ish2  !   "
208     REAL(wp), POINTER, DIMENSION(:)    :: zwork !   " 
209     REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   "
210     !!---------------------------------------------------------------------   
211     !
212     IF( nn_timing == 1 )   CALL timing_start('dia_dct')
213
214     IF( lk_mpp )THEN
215        itotal = nb_sec_max*nb_type_class*nb_class_max
216        CALL wrk_alloc( itotal                                , zwork ) 
217        CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum  )
218     ENDIF   
219 
220     ! Initialise arrays
221     zwork(:) = 0.0 
222     zsum(:,:,:) = 0.0
223
224     IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
225         WRITE(numout,*) " "
226         WRITE(numout,*) "diadct: compute transport"
227         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
228         WRITE(numout,*) "nb_sec = ",nb_sec
229     ENDIF
230
231 
232     ! Compute transport and write only at nn_dctwri
233     IF( MOD(kt,nn_dct)==0 ) THEN
234
235        DO jsec=1,nb_sec
236
237           !debug this section computing ?
238           lldebug=.FALSE.
239           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 ) lldebug=.TRUE. 
240
241           !Compute transport through section 
242           CALL transport(secs(jsec),lldebug,jsec) 
243
244        ENDDO
245             
246        IF( MOD(kt,nn_dctwri)==0 )THEN
247
248           IF( kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: average transports and write at kt = ",kt         
249 
250           !! divide arrays by nn_dctwri/nn_dct to obtain average
251           transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 
252           transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct) 
253 
254           ! Sum over each class
255           DO jsec=1,nb_sec 
256              CALL dia_dct_sum(secs(jsec),jsec) 
257           ENDDO 
258
259           !Sum on all procs
260           IF( lk_mpp )THEN
261              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
262              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
263              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
264              zwork(:)= RESHAPE(zsum(:,:,:), ish )
265              CALL mpp_sum(zwork, ish(1))
266              zsum(:,:,:)= RESHAPE(zwork,ish2)
267              DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
268           ENDIF
269
270           !Write the transport
271           DO jsec=1,nb_sec
272
273              IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
274           
275              !nullify transports values after writing
276              transports_3d(:,jsec,:,:)=0.
277              transports_2d(:,jsec,:  )=0.
278              secs(jsec)%transport(:,:)=0. 
279
280           ENDDO
281
282        ENDIF
283
284     ENDIF
285
286     IF( lk_mpp )THEN
287        itotal = nb_sec_max*nb_type_class*nb_class_max
288        CALL wrk_dealloc( itotal                                , zwork ) 
289        CALL wrk_dealloc( nb_sec_max,nb_type_class,nb_class_max , zsum  )
290     ENDIF   
291
292     IF( nn_timing == 1 )   CALL timing_stop('dia_dct')
293     !
294  END SUBROUTINE dia_dct
295
296  SUBROUTINE readsec 
297     !!---------------------------------------------------------------------
298     !!               ***  ROUTINE readsec  ***
299     !!
300     !!  ** Purpose:
301     !!            Read a binary file(section_ijglobal.diadct)
302     !!            generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
303     !!
304     !!
305     !!---------------------------------------------------------------------
306     !! * Local variables
307     INTEGER :: iptglo , iptloc                               ! Global and local number of points for a section
308     INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2  ! temporary  integer
309     INTEGER :: jsec, jpt                                     ! dummy loop indices
310
311     INTEGER, DIMENSION(2) :: icoord 
312     CHARACTER(len=160)    :: clname                          !filename
313     CHARACTER(len=200)    :: cltmp
314     CHARACTER(len=200)    :: clformat                        !automatic format
315     TYPE(POINT_SECTION),DIMENSION(nb_point_max)  ::coordtemp !contains listpoints coordinates
316                                                              !read in the file
317     INTEGER, POINTER, DIMENSION(:) :: directemp              !contains listpoints directions
318                                                              !read in the files
319     LOGICAL :: llbon                                       ,&!local logical
320                lldebug                                       !debug the section
321     !!-------------------------------------------------------------------------------------
322     CALL wrk_alloc( nb_point_max, directemp )
323
324     !open input file
325     !---------------
326     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
327 
328     !---------------
329     !Read input file
330     !---------------
331     
332     DO jsec=1,nb_sec_max      !loop on the nb_sec sections
333
334        IF (  jsec==nn_secdebug .OR. nn_secdebug==-1  ) &
335           & WRITE(numout,*)'debuging for section number: ',jsec 
336
337        !initialization
338        !---------------
339        secs(jsec)%name=''
340        secs(jsec)%llstrpond    = .FALSE.  ; secs(jsec)%ll_ice_section = .FALSE.
341        secs(jsec)%ll_date_line = .FALSE.  ; secs(jsec)%nb_class       = 0
342        secs(jsec)%zsigi        = 99._wp   ; secs(jsec)%zsigp          = 99._wp
343        secs(jsec)%zsal         = 99._wp   ; secs(jsec)%ztem           = 99._wp
344        secs(jsec)%zlay         = 99._wp         
345        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0
346
347        !read section's number / name / computing choices / classes / slopeSection / points number
348        !-----------------------------------------------------------------------------------------
349        READ(numdct_in,iostat=iost)isec
350        IF (iost .NE. 0 )EXIT !end of file
351        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
352        IF( jsec .NE. isec )  CALL ctl_stop( cltmp )
353
354        IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )WRITE(numout,*)"isec ",isec 
355
356        READ(numdct_in)secs(jsec)%name
357        READ(numdct_in)secs(jsec)%llstrpond
358        READ(numdct_in)secs(jsec)%ll_ice_section
359        READ(numdct_in)secs(jsec)%ll_date_line
360        READ(numdct_in)secs(jsec)%coordSec
361        READ(numdct_in)secs(jsec)%nb_class
362        READ(numdct_in)secs(jsec)%zsigi
363        READ(numdct_in)secs(jsec)%zsigp
364        READ(numdct_in)secs(jsec)%zsal
365        READ(numdct_in)secs(jsec)%ztem
366        READ(numdct_in)secs(jsec)%zlay
367        READ(numdct_in)secs(jsec)%slopeSection
368        READ(numdct_in)iptglo
369
370        !debug
371        !-----
372
373        IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
374         
375            WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))' 
376
377            WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name)
378            WRITE(numout,*)       "      Compute heat and salt transport ? ",secs(jsec)%llstrpond
379            WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section
380            WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line
381            WRITE(numout,*)       "      Slope section :                   ",secs(jsec)%slopeSection
382            WRITE(numout,*)       "      Number of points in the section:  ",iptglo
383            WRITE(numout,*)       "      Number of classes                 ",secs(jsec)%nb_class
384            WRITE(numout,clformat)"      Insitu density classes :          ",secs(jsec)%zsigi
385            WRITE(numout,clformat)"      Potential density classes :       ",secs(jsec)%zsigp
386            WRITE(numout,clformat)"      Salinity classes :                ",secs(jsec)%zsal
387            WRITE(numout,clformat)"      Temperature classes :             ",secs(jsec)%ztem
388            WRITE(numout,clformat)"      Depth classes :                   ",secs(jsec)%zlay
389        ENDIF               
390
391        IF( iptglo /= 0 )THEN
392             
393           !read points'coordinates and directions
394           !--------------------------------------
395           coordtemp(:) = POINT_SECTION(0,0) !list of points read
396           directemp(:) = 0                  !value of directions of each points
397           DO jpt=1,iptglo
398              READ(numdct_in) i1, i2
399              coordtemp(jpt)%I = i1 
400              coordtemp(jpt)%J = i2
401           ENDDO
402           READ(numdct_in) directemp(1:iptglo)
403   
404           !debug
405           !-----
406           IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
407              WRITE(numout,*)"      List of points in global domain:"
408              DO jpt=1,iptglo
409                 WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt),directemp(jpt)
410              ENDDO                 
411           ENDIF
412 
413           !Now each proc selects only points that are in its domain:
414           !--------------------------------------------------------
415           iptloc = 0                    ! initialize number of points selected
416           DO jpt = 1, iptglo            ! loop on listpoint read in the file
417              !     
418              iiglo=coordtemp(jpt)%I          ! global coordinates of the point
419              ijglo=coordtemp(jpt)%J          !  "
420
421              IF( iiglo==jpiglo .AND. nimpp==1 )   iiglo = 2         !!gm BUG: Hard coded periodicity !
422
423              iiloc=iiglo-nimpp+1   ! local coordinates of the point
424              ijloc=ijglo-njmpp+1   !  "
425
426              !verify if the point is on the local domain:(1,nlei)*(1,nlej)
427              IF( iiloc >= 1 .AND. iiloc <= nlei .AND. &
428                  ijloc >= 1 .AND. ijloc <= nlej       )THEN
429                 iptloc = iptloc + 1                                                 ! count local points
430                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
431                 secs(jsec)%direction(iptloc) = directemp(jpt)                       ! store local direction
432              ENDIF
433              !
434           END DO
435     
436           secs(jsec)%nb_point=iptloc !store number of section's points
437
438           !debug
439           !-----
440           IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
441              WRITE(numout,*)"      List of points selected by the proc:"
442              DO jpt = 1,iptloc
443                 iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
444                 ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
445                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
446              ENDDO
447           ENDIF
448
449              IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
450              DO jpt = 1,iptloc
451                 iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
452                 ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
453              ENDDO
454              ENDIF
455
456           !remove redundant points between processors
457           !------------------------------------------
458           lldebug = .FALSE. ; IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) lldebug = .TRUE.
459           IF( iptloc .NE. 0 )THEN
460              CALL removepoints(secs(jsec),'I','top_list',lldebug)
461              CALL removepoints(secs(jsec),'I','bot_list',lldebug)
462              CALL removepoints(secs(jsec),'J','top_list',lldebug)
463              CALL removepoints(secs(jsec),'J','bot_list',lldebug)
464           ENDIF
465           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
466              DO jpt = 1,secs(jsec)%nb_point
467                 iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
468                 ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
469              ENDDO
470           ENDIF
471
472           !debug
473           !-----
474           IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
475              WRITE(numout,*)"      List of points after removepoints:"
476              iptloc = secs(jsec)%nb_point
477              DO jpt = 1,iptloc
478                 iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
479                 ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
480                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
481                 CALL FLUSH(numout)
482              ENDDO
483           ENDIF
484
485        ELSE  ! iptglo = 0
486           IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )&
487              WRITE(numout,*)'   No points for this section.'
488        ENDIF
489
490     ENDDO !end of the loop on jsec
491 
492     nb_sec = jsec-1   !number of section read in the file
493
494     CALL wrk_dealloc( nb_point_max, directemp )
495     !
496  END SUBROUTINE readsec
497
498  SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
499     !!---------------------------------------------------------------------------
500     !!             *** function removepoints
501     !!
502     !!   ** Purpose :: Remove points which are common to 2 procs
503     !!
504     !----------------------------------------------------------------------------
505     !! * arguments
506     TYPE(SECTION),INTENT(INOUT) :: sec
507     CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
508     CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
509     LOGICAL,INTENT(IN)          :: ld_debug                     
510
511     !! * Local variables
512     INTEGER :: iextr         ,& !extremity of listpoint that we verify
513                iind          ,& !coord     of listpoint that we verify
514                itest         ,& !indice value of the side of the domain
515                                 !where points could be redundant
516                isgn          ,& ! isgn= 1 : scan listpoint from start to end
517                                 ! isgn=-1 : scan listpoint from end to start
518                istart,iend      !first and last points selected in listpoint
519     INTEGER :: jpoint           !loop on list points
520     INTEGER, POINTER, DIMENSION(:)   :: idirec !contains temporary sec%direction
521     INTEGER, POINTER, DIMENSION(:,:) :: icoord !contains temporary sec%listpoint
522     !----------------------------------------------------------------------------
523     CALL wrk_alloc(    nb_point_max, idirec )
524     CALL wrk_alloc( 2, nb_point_max, icoord )
525
526     IF( ld_debug )WRITE(numout,*)'      -------------------------'
527     IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
528
529     !iextr=extremity of list_point that we verify
530     IF      ( cdextr=='bot_list' )THEN ; iextr=1            ; isgn=1
531     ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
532     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdextr")
533     ENDIF
534 
535     !which coordinate shall we verify ?
536     IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1
537     ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2
538     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind") 
539     ENDIF
540
541     IF( ld_debug )THEN
542        WRITE(numout,*)'      case: coord/list extr/domain side'
543        WRITE(numout,*)'      ', cdind,' ',cdextr,' ',itest
544        WRITE(numout,*)'      Actual number of points: ',sec%nb_point
545     ENDIF
546
547     icoord(1,1:nb_point_max) = sec%listPoint%I
548     icoord(2,1:nb_point_max) = sec%listPoint%J
549     idirec                   = sec%direction
550     sec%listPoint            = POINT_SECTION(0,0)
551     sec%direction            = 0
552
553     jpoint=iextr+isgn
554     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
555         IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
556         ELSE                                                                               ; EXIT
557         ENDIF
558     ENDDO 
559
560     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
561     ELSE                        ; istart=1        ; iend=jpoint+1
562     ENDIF
563
564     sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
565     sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
566     sec%direction(1:1+iend-istart)   = idirec(istart:iend)
567     sec%nb_point                     = iend-istart+1
568     
569     IF( ld_debug )THEN
570        WRITE(numout,*)'      Number of points after removepoints :',sec%nb_point
571        WRITE(numout,*)'      sec%direction after removepoints :',sec%direction(1:sec%nb_point)
572     ENDIF
573
574     CALL wrk_dealloc(    nb_point_max, idirec )
575     CALL wrk_dealloc( 2, nb_point_max, icoord )
576  END SUBROUTINE removepoints
577
578  SUBROUTINE transport(sec,ld_debug,jsec)
579     !!-------------------------------------------------------------------------------------------
580     !!                     ***  ROUTINE transport  ***
581     !!
582     !!  Purpose ::  Compute the transport for each point in a section
583     !!
584     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
585     !!              Be aware :           
586     !!              One section is a sum of segments
587     !!              One segment is defined by 2 consecutive points in sec%listPoint
588     !!              All points of sec%listPoint are positioned on the F-point of the cell
589     !!
590     !!              There are two loops:                 
591     !!              loop on the segment between 2 nodes
592     !!              loop on the level jk !!
593     !!
594     !!  Output  ::  Arrays containing the volume,density,heat,salt transports for each i
595     !!              point in a section, summed over each nn_dct.
596     !!
597     !!-------------------------------------------------------------------------------------------
598     !! * Arguments
599     TYPE(SECTION),INTENT(INOUT) :: sec
600     LOGICAL      ,INTENT(IN)    :: ld_debug
601     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
602   
603     !! * Local variables
604     INTEGER             :: jk, jseg, jclass,jl,                 &!loop on level/segment/classes/ice categories
605                            isgnu, isgnv                          !
606     REAL(wp)            :: zumid, zvmid,                        &!U/V velocity on a cell segment
607                            zumid_ice, zvmid_ice,                &!U/V ice velocity
608                            zTnorm                                !transport of velocity through one cell's sides
609     REAL(wp)            :: ztn, zsn, zrhoi, zrhop, zsshn, zdep !temperature/salinity/potential density/ssh/depth at u/v point
610
611     TYPE(POINT_SECTION) :: k
612      !!--------------------------------------------------------
613      !
614      IF( ld_debug )WRITE(numout,*)'      Compute transport'
615
616      !---------------------------!
617      !  COMPUTE TRANSPORT        !
618      !---------------------------!
619      IF(sec%nb_point .NE. 0)THEN   
620
621         !----------------------------------------------------------------------------------------------------
622         !Compute sign for velocities:
623         !
624         !convention:
625         !   non horizontal section: direction + is toward left hand of section
626         !       horizontal section: direction + is toward north of section
627         !
628         !
629         !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
630         !       ----------------      -----------------     ---------------             --------------
631         !
632         !   isgnv=1         direction +     
633         !  ______         _____             ______                                                   
634         !        |           //|            |                  |                         direction +   
635         !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
636         !        |_______  //         ______|    \\            | ---\                        |
637         !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
638         !               |             |          __\\|         |                   
639         !               |             |     direction +        |                      isgnv=1                                 
640         !                                                     
641         !----------------------------------------------------------------------------------------------------
642         isgnu = 1
643         IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
644         ELSE                                ; isgnv =  1
645         ENDIF
646         IF( sec%slopeSection .GE. 9999. )     isgnv =  1
647
648         IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
649
650         !--------------------------------------!
651         ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
652         !--------------------------------------!
653         DO jseg=1,MAX(sec%nb_point-1,0)
654             
655            !-------------------------------------------------------------------------------------------
656            ! Select the appropriate coordinate for computing the velocity of the segment
657            !
658            !                      CASE(0)                                    Case (2)
659            !                      -------                                    --------
660            !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
661            !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
662            !                                                                            |
663            !                                                                            |
664            !                                                                            |
665            !                      Case (3)                                            U(i,j)
666            !                      --------                                              |
667            !                                                                            |
668            !  listPoint(jseg+1) F(i,j+1)                                                |
669            !                        |                                                   |
670            !                        |                                                   |
671            !                        |                                 listPoint(jseg+1) F(i,j-1)
672            !                        |                                           
673            !                        |                                           
674            !                     U(i,j+1)                                           
675            !                        |                                       Case(1)     
676            !                        |                                       ------     
677            !                        |                                           
678            !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
679            !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
680            ! listPoint(jseg)     F(i,j)
681            !
682            !-------------------------------------------------------------------------------------------
683
684            SELECT CASE( sec%direction(jseg) )
685            CASE(0)   ;    k = sec%listPoint(jseg)
686            CASE(1)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
687            CASE(2)   ;    k = sec%listPoint(jseg)
688            CASE(3)   ;    k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
689            END SELECT
690
691            !---------------------------|
692            !     LOOP ON THE LEVEL     |
693            !---------------------------|
694            DO jk = 1, mbkt(k%I,k%J)            !Sum of the transport on the vertical
695            !           ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
696            SELECT CASE( sec%direction(jseg) )
697               CASE(0,1) 
698                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
699                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
700                  zrhop = interp(k%I,k%J,jk,'V',rhop) 
701                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
702                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
703               CASE(2,3) 
704                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
705                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
706                  zrhop = interp(k%I,k%J,jk,'U',rhop) 
707                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
708                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
709               END SELECT 
710               !
711               zdep= gdept_n(k%I,k%J,jk) 
712 
713               SELECT CASE( sec%direction(jseg) )                !compute velocity with the correct direction
714               CASE(0,1)   
715                  zumid=0._wp
716                  zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
717               CASE(2,3) 
718                  zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
719                  zvmid=0._wp
720               END SELECT 
721 
722               !zTnorm=transport through one cell;
723               !velocity* cell's length * cell's thickness
724               zTnorm = zumid*e2u(k%I,k%J) * e3u_n(k%I,k%J,jk)     & 
725                  &   + zvmid*e1v(k%I,k%J) * e3v_n(k%I,k%J,jk) 
726
727!!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!!
728               IF( ln_linssh ) THEN              !add transport due to free surface
729                  IF( jk==1 ) THEN
730                     zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk)   & 
731                        &            + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
732                  ENDIF
733               ENDIF
734!!gm end
735              !COMPUTE TRANSPORT 
736 
737              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
738 
739              IF( sec%llstrpond ) THEN
740                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp
741                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001
742              ENDIF
743   
744           END DO !end of loop on the level
745
746#if defined key_lim3
747
748           !ICE CASE   
749           !------------
750           IF( sec%ll_ice_section )THEN
751              SELECT CASE (sec%direction(jseg))
752              CASE(0)
753                 zumid_ice = 0
754                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
755              CASE(1)
756                 zumid_ice = 0
757                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
758              CASE(2)
759                 zvmid_ice = 0
760                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
761              CASE(3)
762                 zvmid_ice = 0
763                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
764              END SELECT
765   
766              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
767
768#if defined key_lim3
769              DO jl=1,jpl
770                 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*       &
771                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) *  &
772                                  ( h_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  &
773                                    h_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
774                                   
775                 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
776                                    a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
777              END DO
778#endif
779   
780           ENDIF !end of ice case
781#endif
782 
783        END DO !end of loop on the segment
784
785     ENDIF !end of sec%nb_point =0 case
786     !
787  END SUBROUTINE transport
788
789
790  SUBROUTINE dia_dct_sum(sec,jsec) 
791     !!-------------------------------------------------------------
792     !! Purpose: Average the transport over nn_dctwri time steps 
793     !! and sum over the density/salinity/temperature/depth classes
794     !!
795     !! Method:   Sum over relevant grid cells to obtain values 
796     !!           for each class
797     !!              There are several loops:                 
798     !!              loop on the segment between 2 nodes
799     !!              loop on the level jk
800     !!              loop on the density/temperature/salinity/level classes
801     !!              test on the density/temperature/salinity/level
802     !!
803     !!  Note:    Transport through a given section is equal to the sum of transports
804     !!           computed on each proc.
805     !!           On each proc,transport is equal to the sum of transport computed through
806     !!           segments linking each point of sec%listPoint  with the next one.   
807     !!
808     !!-------------------------------------------------------------
809     TYPE(SECTION),INTENT(INOUT) :: sec 
810     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
811 
812     TYPE(POINT_SECTION) :: k 
813     INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes 
814     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point
815     !!-------------------------------------------------------------
816 
817     !! Sum the relevant segments to obtain values for each class
818     IF(sec%nb_point .NE. 0)THEN   
819 
820        !--------------------------------------!
821        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
822        !--------------------------------------!
823        DO jseg=1,MAX(sec%nb_point-1,0) 
824           
825           !-------------------------------------------------------------------------------------------
826           ! Select the appropriate coordinate for computing the velocity of the segment
827           !
828           !                      CASE(0)                                    Case (2)
829           !                      -------                                    --------
830           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
831           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
832           !                                                                            |
833           !                                                                            |
834           !                                                                            |
835           !                      Case (3)                                            U(i,j)
836           !                      --------                                              |
837           !                                                                            |
838           !  listPoint(jseg+1) F(i,j+1)                                                |
839           !                        |                                                   |
840           !                        |                                                   |
841           !                        |                                 listPoint(jseg+1) F(i,j-1)
842           !                        |                                             
843           !                        |                                             
844           !                     U(i,j+1)                                             
845           !                        |                                       Case(1)     
846           !                        |                                       ------       
847           !                        |                                             
848           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
849           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
850           ! listPoint(jseg)     F(i,j)
851           
852           !-------------------------------------------------------------------------------------------
853 
854           SELECT CASE( sec%direction(jseg) ) 
855           CASE(0)  ;   k = sec%listPoint(jseg) 
856           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
857           CASE(2)  ;   k = sec%listPoint(jseg) 
858           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
859           END SELECT 
860 
861           !---------------------------|
862           !     LOOP ON THE LEVEL     |
863           !---------------------------|
864           !Sum of the transport on the vertical 
865           DO jk=1,mbkt(k%I,k%J) 
866 
867              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
868              SELECT CASE( sec%direction(jseg) ) 
869              CASE(0,1) 
870                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
871                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
872                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
873                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
874
875              CASE(2,3) 
876                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
877                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
878                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
879                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
880                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
881              END SELECT
882 
883              zdep= gdept_n(k%I,k%J,jk) 
884 
885              !-------------------------------
886              !  LOOP ON THE DENSITY CLASSES |
887              !-------------------------------
888              !The computation is made for each density/temperature/salinity/depth class
889              DO jclass=1,MAX(1,sec%nb_class-1) 
890 
891                 !----------------------------------------------!
892                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 
893                 !----------------------------------------------!
894
895                 IF ( (                                                    & 
896                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
897                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
898                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
899 
900                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
901                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
902                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
903 
904                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
905                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
906                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
907 
908                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
909                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
910                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
911 
912                    ((( zdep .GE. sec%zlay(jclass)) .AND.                & 
913                    (   zdep .LE. sec%zlay(jclass+1))) .OR.              & 
914                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
915                                                                   ))   THEN 
916 
917                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
918                    !----------------------------------------------------------------------------
919                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN 
920                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
921                    ELSE
922                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
923                    ENDIF
924                    IF( sec%llstrpond )THEN
925 
926                       IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN
927                          sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk) 
928                       ELSE
929                          sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk) 
930                       ENDIF
931 
932                       IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN
933                          sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk) 
934                       ELSE
935                          sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk) 
936                       ENDIF
937 
938                    ELSE
939                       sec%transport( 3,jclass) = 0._wp 
940                       sec%transport( 4,jclass) = 0._wp 
941                       sec%transport( 5,jclass) = 0._wp 
942                       sec%transport( 6,jclass) = 0._wp 
943                    ENDIF
944 
945                 ENDIF ! end of test if point is in class
946   
947              END DO ! end of loop on the classes
948 
949           END DO ! loop over jk
950 
951#if defined key_lim3 
952 
953           !ICE CASE     
954           IF( sec%ll_ice_section )THEN
955 
956              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
957                 sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6 
958              ELSE
959                 sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6 
960              ENDIF
961 
962              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
963                 sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6 
964              ELSE
965                 sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6 
966              ENDIF
967 
968           ENDIF !end of ice case
969#endif
970 
971        END DO !end of loop on the segment
972 
973     ELSE  !if sec%nb_point =0
974        sec%transport(1:2,:)=0. 
975        IF (sec%llstrpond) sec%transport(3:6,:)=0. 
976        IF (sec%ll_ice_section) sec%transport(7:10,:)=0. 
977     ENDIF !end of sec%nb_point =0 case
978 
979  END SUBROUTINE dia_dct_sum 
980
981
982  SUBROUTINE dia_dct_wri(kt,ksec,sec)
983     !!-------------------------------------------------------------
984     !! Write transport output in numdct
985     !!
986     !! Purpose: Write  transports in ascii files
987     !!
988     !! Method:
989     !!        1. Write volume transports in "volume_transport"
990     !!           Unit: Sv : area * Velocity / 1.e6
991     !!
992     !!        2. Write heat transports in "heat_transport"
993     !!           Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
994     !!
995     !!        3. Write salt transports in "salt_transport"
996     !!           Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9
997     !!
998     !!-------------------------------------------------------------
999     !!arguments
1000     INTEGER, INTENT(IN)          :: kt          ! time-step
1001     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
1002     INTEGER ,INTENT(IN)          :: ksec        ! section number
1003
1004     !!local declarations
1005     INTEGER               :: jclass             ! Dummy loop
1006     CHARACTER(len=2)      :: classe             ! Classname
1007     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
1008     REAL(wp)              :: zslope             ! section's slope coeff
1009     !
1010     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
1011     !!-------------------------------------------------------------
1012     CALL wrk_alloc(nb_type_class , zsumclasses ) 
1013
1014     zsumclasses(:)=0._wp
1015     zslope = sec%slopeSection       
1016
1017 
1018     DO jclass=1,MAX(1,sec%nb_class-1)
1019
1020        classe   = 'N       '
1021        zbnd1   = 0._wp
1022        zbnd2   = 0._wp
1023        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
1024
1025   
1026        !insitu density classes transports
1027        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
1028            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
1029           classe = 'DI       '
1030           zbnd1 = sec%zsigi(jclass)
1031           zbnd2 = sec%zsigi(jclass+1)
1032        ENDIF
1033        !potential density classes transports
1034        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
1035            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
1036           classe = 'DP      '
1037           zbnd1 = sec%zsigp(jclass)
1038           zbnd2 = sec%zsigp(jclass+1)
1039        ENDIF
1040        !depth classes transports
1041        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
1042            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
1043           classe = 'Z       '
1044           zbnd1 = sec%zlay(jclass)
1045           zbnd2 = sec%zlay(jclass+1)
1046        ENDIF
1047        !salinity classes transports
1048        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
1049            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
1050           classe = 'S       '
1051           zbnd1 = sec%zsal(jclass)
1052           zbnd2 = sec%zsal(jclass+1)   
1053        ENDIF
1054        !temperature classes transports
1055        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
1056            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
1057           classe = 'T       '
1058           zbnd1 = sec%ztem(jclass)
1059           zbnd2 = sec%ztem(jclass+1)
1060        ENDIF
1061                 
1062        !write volume transport per class
1063        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1064                              jclass,classe,zbnd1,zbnd2,&
1065                              sec%transport(1,jclass),sec%transport(2,jclass), &
1066                              sec%transport(1,jclass)+sec%transport(2,jclass)
1067
1068        IF( sec%llstrpond )THEN
1069
1070           !write heat transport per class:
1071           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  &
1072                              jclass,classe,zbnd1,zbnd2,&
1073                              sec%transport(3,jclass)*1.e-15,sec%transport(4,jclass)*1.e-15, &
1074                              ( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15
1075           !write salt transport per class
1076           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  &
1077                              jclass,classe,zbnd1,zbnd2,&
1078                              sec%transport(5,jclass)*1.e-9,sec%transport(6,jclass)*1.e-9,&
1079                              (sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9
1080        ENDIF
1081
1082     ENDDO
1083
1084     zbnd1 = 0._wp
1085     zbnd2 = 0._wp
1086     jclass=0
1087
1088     !write total volume transport
1089     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1090                           jclass,"total",zbnd1,zbnd2,&
1091                           zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
1092
1093     IF( sec%llstrpond )THEN
1094
1095        !write total heat transport
1096        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1097                           jclass,"total",zbnd1,zbnd2,&
1098                           zsumclasses(3)*1.e-15,zsumclasses(4)*1.e-15,&
1099                           (zsumclasses(3)+zsumclasses(4) )*1.e-15
1100        !write total salt transport
1101        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1102                           jclass,"total",zbnd1,zbnd2,&
1103                           zsumclasses(5)*1.e-9,zsumclasses(6)*1.e-9,&
1104                           (zsumclasses(5)+zsumclasses(6))*1.e-9
1105     ENDIF
1106
1107     
1108     IF ( sec%ll_ice_section) THEN
1109        !write total ice volume transport
1110        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1111                              jclass,"ice_vol",zbnd1,zbnd2,&
1112                              sec%transport(7,1),sec%transport(8,1),&
1113                              sec%transport(7,1)+sec%transport(8,1)
1114        !write total ice surface transport
1115        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1116                              jclass,"ice_surf",zbnd1,zbnd2,&
1117                              sec%transport(9,1),sec%transport(10,1), &
1118                              sec%transport(9,1)+sec%transport(10,1) 
1119      ENDIF
1120                                             
1121118   FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
1122119   FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
1123
1124      CALL wrk_dealloc(nb_type_class , zsumclasses ) 
1125      !
1126   END SUBROUTINE dia_dct_wri
1127
1128
1129   FUNCTION interp(ki, kj, kk, cd_point, ptab)
1130  !!----------------------------------------------------------------------
1131  !!
1132  !!   Purpose: compute temperature/salinity/density at U-point or V-point
1133  !!   --------
1134  !!
1135  !!   Method:
1136  !!   ------
1137  !!
1138  !!   ====> full step and partial step
1139  !!
1140  !!
1141  !!    |    I          |    I+1           |    Z=temperature/salinity/density at U-poinT
1142  !!    |               |                  |
1143  !!  ----------------------------------------  1. Veritcal interpolation: compute zbis
1144  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
1145  !!    |               |                  |       zbis =
1146  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
1147  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
1148  !!    |               |                  |
1149  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
1150  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K) 
1151  !!    |     .         |                  |
1152  !!    |     .         |                  |       interp = ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)
1153  !!    |     .         |                  |
1154  !!  ------------------------------------------
1155  !!    |     .         |                  |
1156  !!    |     .         |                  |
1157  !!    |     .         |                  |
1158  !!K   |    zbis.......U...ptab(I+1,J,K)  |
1159  !!    |     .         |                  |
1160  !!    | ptab(I,J,K)   |                  |
1161  !!    |               |------------------|
1162  !!    |               | partials         |
1163  !!    |               |  steps           |
1164  !!  -------------------------------------------
1165  !!    <----zet1------><----zet2--------->
1166  !!
1167  !!
1168  !!   ====>  s-coordinate
1169  !!     
1170  !!    |                |                  |   1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
1171  !!    |                |                  |      Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
1172  !!    |                | ptab(I+1,J,K)    |
1173  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point
1174  !!    |                |      ^           |   
1175  !!    |                |      | zdep2     |   
1176  !!    |                |      |           |   
1177  !!    |       ^        U      v           |
1178  !!    |       |        |                  |
1179  !!    |       | zdep1  |                  |   
1180  !!    |       v        |                  |
1181  !!    |      T1        |                  |
1182  !!    | ptab(I,J,K)    |                  |
1183  !!    |                |                  |
1184  !!    |                |                  |
1185  !!
1186  !!    <----zet1--------><----zet2--------->
1187  !!
1188  !!----------------------------------------------------------------------
1189  !*arguments
1190  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
1191  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
1192  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk )
1193  REAL(wp)                                     :: interp       ! interpolated variable
1194
1195  !*local declations
1196  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
1197  REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu            ! local real
1198  REAL(wp):: zet1, zet2                                        ! weight for interpolation
1199  REAL(wp):: zdep1,zdep2                                       ! differences of depth
1200  REAL(wp):: zmsk                                              ! mask value
1201  !!----------------------------------------------------------------------
1202
1203  IF( cd_point=='U' )THEN
1204     ii1 = ki    ; ij1 = kj 
1205     ii2 = ki+1  ; ij2 = kj 
1206
1207     zet1=e1t(ii1,ij1)
1208     zet2=e1t(ii2,ij2)
1209     zmsk=umask(ii1,ij1,kk)
1210 
1211
1212  ELSE ! cd_point=='V'
1213     ii1 = ki    ; ij1 = kj 
1214     ii2 = ki    ; ij2 = kj+1 
1215
1216     zet1=e2t(ii1,ij1)
1217     zet2=e2t(ii2,ij2)
1218     zmsk=vmask(ii1,ij1,kk)
1219
1220  ENDIF
1221
1222  IF( ln_sco )THEN   ! s-coordinate case
1223
1224     zdepu = ( gdept_n(ii1,ij1,kk) +  gdept_n(ii2,ij2,kk) ) * 0.5_wp 
1225     zdep1 = gdept_n(ii1,ij1,kk) - zdepu
1226     zdep2 = gdept_n(ii2,ij2,kk) - zdepu
1227
1228     ! weights
1229     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1230     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1231 
1232     ! result
1233     interp = zmsk * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )   
1234
1235
1236  ELSE       ! full step or partial step case
1237
1238     ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk) 
1239     zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
1240     zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
1241
1242     IF(kk .NE. 1)THEN
1243
1244        IF( ze3t >= 0. )THEN 
1245           ! zbis
1246           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
1247           ! result
1248            interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
1249        ELSE
1250           ! zbis
1251           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1252           ! result
1253           interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1254        ENDIF   
1255
1256     ELSE
1257        interp = zmsk * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1258     ENDIF
1259
1260  ENDIF
1261      !
1262   END FUNCTION interp
1263
1264#else
1265   !!----------------------------------------------------------------------
1266   !!   Default option :                                       Dummy module
1267   !!----------------------------------------------------------------------
1268   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
1269   PUBLIC 
1270   !! $Id$
1271CONTAINS
1272
1273   SUBROUTINE dia_dct_init          ! Dummy routine
1274      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
1275   END SUBROUTINE dia_dct_init
1276
1277   SUBROUTINE dia_dct( kt )         ! Dummy routine
1278      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
1279      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
1280   END SUBROUTINE dia_dct
1281#endif
1282
1283   !!======================================================================
1284END MODULE diadct
Note: See TracBrowser for help on using the repository browser.