New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diadct.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diadct.F90 @ 12353

Last change on this file since 12353 was 11960, checked in by acc, 5 years ago
Branch 2019/dev_r11943_MERGE_2019. Merge in changes from 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. (svn merge -r 11614:11954). Resolved tree conflicts and one actual conflict. Sette tested(these changes alter the ext/AGRIF reference; remember to update). See ticket #2341
Property svn:keywords set to `Id`
File size: 58.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: