source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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