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

Last change on this file since 9124 was 9124, checked in by gm, 4 years ago

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

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