source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 5726

Last change on this file since 5726 was 5726, checked in by jpalmier, 5 years ago

JPALM — 10-09-2015 — add MEDUSA in the branch ; adapted TOP_SRC to MEDUSA ; remove some svn keywords in the branch

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