Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

diadct.F90 @ 9570

Last change on this file since 9570 was 9570, checked in by nicolasmartin, 6 years ago

Global renaming for core routines (./NEMO)

Folders
- LIM_SRC_3 -> ICE_SRC
- OPA_SRC -> OCE_SRC
CPP key: key_lim3 -> key_si3
Modules, (sub)routines and variables names
- MPI: mpi_comm_opa -> mpi_comm_oce, MPI_COMM_OPA -> MPI_COMM_OCE, mpi_init_opa -> mpi_init_oce
- AGRIF: agrif_opa_* -> agrif_oce_*, agrif_lim3_* -> agrif_si3_* and few more
- TOP-PISCES: p.zlim -> p.zice, namp.zlim -> namp.zice
Comments
- NEMO/OPA -> NEMO/OCE
- ESIM|LIM3 -> SI3

Property svn:keywords set to Id

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

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

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

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OCE_SRC/DIA/diadct.F90 @ 9570

Download in other formats: