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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diadct.F90 @ 12067

Last change on this file since 12067 was 11822, checked in by acc, 5 years ago
Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.
Property svn:keywords set to `Id`
File size: 58.6 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	!!----------------------------------------------------------------------
[11822]	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
[11822]	42	PUBLIC dia_dct_init ! routine called by nemogcm.F90
[2848]	43
[11822]	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
[11822]	99	INTEGER FUNCTION diadct_alloc()
	100	!!----------------------------------------------------------------------
	101	!! * FUNCTION diadct_alloc *
	102	!!----------------------------------------------------------------------
[3680]	103
[11822]	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
[11822]	107	CALL mpp_sum( 'diadct', diadct_alloc )
	108	IF( diadct_alloc /= 0 ) CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )
[3680]	109
[11822]	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	!!
[11822]	122	NAMELIST/nam_diadct/ln_diadct, nn_dct, nn_dctwri, nn_secdebug
[9124]	123	!!---------------------------------------------------------------------
[2848]	124
[11822]	125	REWIND( numnam_ref ) ! Namelist nam_diadct in reference namelist : Diagnostic: transport through sections
	126	READ ( numnam_ref, nam_diadct, IOSTAT = ios, ERR = 901)
	127	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadct in reference namelist' )
[2848]	128
[11822]	129	REWIND( numnam_cfg ) ! Namelist nam_diadct in configuration namelist : Diagnostic: transport through sections
	130	READ ( numnam_cfg, nam_diadct, IOSTAT = ios, ERR = 902 )
	131	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_diadct in configuration namelist' )
	132	IF(lwm) WRITE ( numond, nam_diadct )
[4147]	133
[2848]	134	IF( lwp ) THEN
	135	WRITE(numout,*) " "
	136	WRITE(numout,*) "diadct_init: compute transports through sections "
	137	WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
[11822]	138	WRITE(numout,*) " Calculate transport thru sections: ln_diadct = ", ln_diadct
	139	WRITE(numout,*) " Frequency of computation: nn_dct = ", nn_dct
	140	WRITE(numout,*) " Frequency of write: nn_dctwri = ", nn_dctwri
[2848]	141
	142	IF ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
	143	WRITE(numout,*)" Debug section number: ", nn_secdebug
	144	ELSE IF ( nn_secdebug == 0 )THEN ; WRITE(numout,*)" No section to debug"
	145	ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)" Debug all sections"
	146	ELSE ; WRITE(numout,*)" Wrong value for nn_secdebug : ",nn_secdebug
	147	ENDIF
[11822]	148	ENDIF
[2848]	149
[11822]	150	IF( ln_diadct ) THEN
	151	! control
[2848]	152	IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0) &
[11822]	153	& CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
[2848]	154
[11822]	155	! allocate dia_dct arrays
	156	IF( diadct_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )
[2848]	157
[11822]	158	!Read section_ijglobal.diadct
	159	CALL readsec
[2848]	160
[11822]	161	!open output file
	162	IF( lwm ) THEN
	163	CALL ctl_opn( numdct_vol, 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
	164	CALL ctl_opn( numdct_heat, 'heat_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
	165	CALL ctl_opn( numdct_salt, 'salt_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
	166	ENDIF
	167
	168	! Initialise arrays to zero
	169	transports_3d(:,:,:,:)=0.0
	170	transports_2d(:,:,:) =0.0
	171	!
[2927]	172	ENDIF
[3294]	173	!
[2848]	174	END SUBROUTINE dia_dct_init
	175
	176
[10965]	177	SUBROUTINE dia_dct( kt, Kmm )
[2848]	178	!!---------------------------------------------------------------------
	179	!! * ROUTINE diadct *
	180	!!
[3680]	181	!! Purpose :: Compute section transports and write it in numdct files
	182	!!
	183	!! Method :: All arrays initialised to zero in dct_init
	184	!! Each nn_dct time step call subroutine 'transports' for
	185	!! each section to sum the transports over each grid cell.
	186	!! Each nn_dctwri time step:
	187	!! Divide the arrays by the number of summations to gain
	188	!! an average value
	189	!! Call dia_dct_sum to sum relevant grid boxes to obtain
	190	!! totals for each class (density, depth, temp or sal)
	191	!! Call dia_dct_wri to write the transports into file
	192	!! Reinitialise all relevant arrays to zero
[2848]	193	!!---------------------------------------------------------------------
[10965]	194	INTEGER, INTENT(in) :: kt ! ocean time step
	195	INTEGER, INTENT(in) :: Kmm ! time level index
[6140]	196	!
[9124]	197	INTEGER :: jsec ! loop on sections
	198	INTEGER :: itotal ! nb_sec_maxnb_type_classnb_class_max
	199	LOGICAL :: lldebug =.FALSE. ! debug a section
	200	INTEGER , DIMENSION(1) :: ish ! work array for mpp_sum
	201	INTEGER , DIMENSION(3) :: ish2 ! "
	202	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zwork ! "
	203	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:):: zsum ! "
[2848]	204	!!---------------------------------------------------------------------
[6140]	205	!
[9124]	206	IF( ln_timing ) CALL timing_start('dia_dct')
[2848]	207
[3294]	208	IF( lk_mpp )THEN
	209	itotal = nb_sec_maxnb_type_classnb_class_max
[9124]	210	ALLOCATE( zwork(itotal) , zsum(nb_sec_max,nb_type_class,nb_class_max) )
[3294]	211	ENDIF
	212
[3680]	213	! Initialise arrays
	214	zwork(:) = 0.0
	215	zsum(:,:,:) = 0.0
	216
[2848]	217	IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
	218	WRITE(numout,*) " "
	219	WRITE(numout,*) "diadct: compute transport"
	220	WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
	221	WRITE(numout,*) "nb_sec = ",nb_sec
	222	ENDIF
	223
	224
	225	! Compute transport and write only at nn_dctwri
	226	IF( MOD(kt,nn_dct)==0 ) THEN
	227
	228	DO jsec=1,nb_sec
	229
	230	!debug this section computing ?
	231	lldebug=.FALSE.
[8126]	232	IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct-1 ) lldebug=.TRUE.
[2848]	233
	234	!Compute transport through section
[10965]	235	CALL transport(Kmm,secs(jsec),lldebug,jsec)
[2908]	236
	237	ENDDO
[2848]	238
[2908]	239	IF( MOD(kt,nn_dctwri)==0 )THEN
[2848]	240
[8126]	241	IF( kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: average transports and write at kt = ",kt
[2848]	242
[3680]	243	!! divide arrays by nn_dctwri/nn_dct to obtain average
	244	transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct)
	245	transports_2d(:,:,:) =transports_2d(:,:,:) /(nn_dctwri/nn_dct)
	246
	247	! Sum over each class
	248	DO jsec=1,nb_sec
[10965]	249	CALL dia_dct_sum(Kmm,secs(jsec),jsec)
[3680]	250	ENDDO
	251
[2908]	252	!Sum on all procs
	253	IF( lk_mpp )THEN
	254	ish(1) = nb_sec_maxnb_type_classnb_class_max
	255	ish2 = (/nb_sec_max,nb_type_class,nb_class_max/)
	256	DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
	257	zwork(:)= RESHAPE(zsum(:,:,:), ish )
[10425]	258	CALL mpp_sum('diadct', zwork, ish(1))
[2908]	259	zsum(:,:,:)= RESHAPE(zwork,ish2)
	260	DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
	261	ENDIF
	262
	263	!Write the transport
	264	DO jsec=1,nb_sec
	265
[5505]	266	IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
[2848]	267
	268	!nullify transports values after writing
[3680]	269	transports_3d(:,jsec,:,:)=0.
	270	transports_2d(:,jsec,: )=0.
[2848]	271	secs(jsec)%transport(:,:)=0.
	272
[2908]	273	ENDDO
	274
	275	ENDIF
	276
[2848]	277	ENDIF
	278
[3294]	279	IF( lk_mpp )THEN
	280	itotal = nb_sec_maxnb_type_classnb_class_max
[9124]	281	DEALLOCATE( zwork , zsum )
[3294]	282	ENDIF
	283
[9124]	284	IF( ln_timing ) CALL timing_stop('dia_dct')
[3294]	285	!
[2848]	286	END SUBROUTINE dia_dct
	287
[9124]	288
[2848]	289	SUBROUTINE readsec
	290	!!---------------------------------------------------------------------
	291	!! * ROUTINE readsec *
	292	!!
[2854]	293	!! ** Purpose:
	294	!! Read a binary file(section_ijglobal.diadct)
	295	!! generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
	296	!!
	297	!!
[2848]	298	!!---------------------------------------------------------------------
	299	INTEGER :: iptglo , iptloc ! Global and local number of points for a section
	300	INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2 ! temporary integer
	301	INTEGER :: jsec, jpt ! dummy loop indices
	302	INTEGER, DIMENSION(2) :: icoord
[9124]	303	LOGICAL :: llbon, lldebug ! local logical
	304	CHARACTER(len=160) :: clname ! filename
[2927]	305	CHARACTER(len=200) :: cltmp
[9124]	306	CHARACTER(len=200) :: clformat !automatic format
	307	TYPE(POINT_SECTION),DIMENSION(nb_point_max) ::coordtemp !contains listpoints coordinates read in the file
	308	INTEGER, DIMENSION(nb_point_max) :: directemp !contains listpoints directions read in the files
[2848]	309	!!-------------------------------------------------------------------------------------
	310
	311	!open input file
	312	!---------------
	313	CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
	314
	315	!---------------
	316	!Read input file
	317	!---------------
	318
	319	DO jsec=1,nb_sec_max !loop on the nb_sec sections
	320
[8126]	321	IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) &
[2848]	322	& WRITE(numout,*)'debuging for section number: ',jsec
	323
	324	!initialization
	325	!---------------
	326	secs(jsec)%name=''
	327	secs(jsec)%llstrpond = .FALSE. ; secs(jsec)%ll_ice_section = .FALSE.
	328	secs(jsec)%ll_date_line = .FALSE. ; secs(jsec)%nb_class = 0
	329	secs(jsec)%zsigi = 99._wp ; secs(jsec)%zsigp = 99._wp
	330	secs(jsec)%zsal = 99._wp ; secs(jsec)%ztem = 99._wp
	331	secs(jsec)%zlay = 99._wp
	332	secs(jsec)%transport = 0._wp ; secs(jsec)%nb_point = 0
	333
	334	!read section's number / name / computing choices / classes / slopeSection / points number
	335	!-----------------------------------------------------------------------------------------
	336	READ(numdct_in,iostat=iost)isec
	337	IF (iost .NE. 0 )EXIT !end of file
[2912]	338	WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
[2848]	339	IF( jsec .NE. isec ) CALL ctl_stop( cltmp )
	340
[8126]	341	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )WRITE(numout,*)"isec ",isec
[2848]	342
	343	READ(numdct_in)secs(jsec)%name
	344	READ(numdct_in)secs(jsec)%llstrpond
	345	READ(numdct_in)secs(jsec)%ll_ice_section
	346	READ(numdct_in)secs(jsec)%ll_date_line
	347	READ(numdct_in)secs(jsec)%coordSec
	348	READ(numdct_in)secs(jsec)%nb_class
	349	READ(numdct_in)secs(jsec)%zsigi
	350	READ(numdct_in)secs(jsec)%zsigp
	351	READ(numdct_in)secs(jsec)%zsal
	352	READ(numdct_in)secs(jsec)%ztem
	353	READ(numdct_in)secs(jsec)%zlay
	354	READ(numdct_in)secs(jsec)%slopeSection
	355	READ(numdct_in)iptglo
	356
	357	!debug
	358	!-----
[2927]	359
[8126]	360	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
[2927]	361
	362	WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))'
	363
	364	WRITE(numout,*) " Section name : ",TRIM(secs(jsec)%name)
	365	WRITE(numout,*) " Compute heat and salt transport ? ",secs(jsec)%llstrpond
	366	WRITE(numout,*) " Compute ice transport ? ",secs(jsec)%ll_ice_section
	367	WRITE(numout,*) " Section crosses date-line ? ",secs(jsec)%ll_date_line
	368	WRITE(numout,*) " Slope section : ",secs(jsec)%slopeSection
	369	WRITE(numout,*) " Number of points in the section: ",iptglo
	370	WRITE(numout,*) " Number of classes ",secs(jsec)%nb_class
	371	WRITE(numout,clformat)" Insitu density classes : ",secs(jsec)%zsigi
	372	WRITE(numout,clformat)" Potential density classes : ",secs(jsec)%zsigp
	373	WRITE(numout,clformat)" Salinity classes : ",secs(jsec)%zsal
	374	WRITE(numout,clformat)" Temperature classes : ",secs(jsec)%ztem
	375	WRITE(numout,clformat)" Depth classes : ",secs(jsec)%zlay
[2848]	376	ENDIF
	377
[7646]	378	IF( iptglo /= 0 )THEN
[2848]	379
	380	!read points'coordinates and directions
	381	!--------------------------------------
	382	coordtemp(:) = POINT_SECTION(0,0) !list of points read
	383	directemp(:) = 0 !value of directions of each points
	384	DO jpt=1,iptglo
[7646]	385	READ(numdct_in) i1, i2
[2848]	386	coordtemp(jpt)%I = i1
	387	coordtemp(jpt)%J = i2
	388	ENDDO
[7646]	389	READ(numdct_in) directemp(1:iptglo)
[2848]	390
	391	!debug
	392	!-----
[8126]	393	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
[2848]	394	WRITE(numout,*)" List of points in global domain:"
	395	DO jpt=1,iptglo
[3632]	396	WRITE(numout,*)' # I J ',jpt,coordtemp(jpt),directemp(jpt)
[2848]	397	ENDDO
	398	ENDIF
	399
	400	!Now each proc selects only points that are in its domain:
	401	!--------------------------------------------------------
[7646]	402	iptloc = 0 ! initialize number of points selected
	403	DO jpt = 1, iptglo ! loop on listpoint read in the file
	404	!
[2848]	405	iiglo=coordtemp(jpt)%I ! global coordinates of the point
	406	ijglo=coordtemp(jpt)%J ! "
	407
[7646]	408	IF( iiglo==jpiglo .AND. nimpp==1 ) iiglo = 2 !!gm BUG: Hard coded periodicity !
[2848]	409
[7646]	410	iiloc=iiglo-nimpp+1 ! local coordinates of the point
	411	ijloc=ijglo-njmpp+1 ! "
[2848]	412
	413	!verify if the point is on the local domain:(1,nlei)*(1,nlej)
[7646]	414	IF( iiloc >= 1 .AND. iiloc <= nlei .AND. &
	415	ijloc >= 1 .AND. ijloc <= nlej )THEN
[2848]	416	iptloc = iptloc + 1 ! count local points
	417	secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
[2912]	418	secs(jsec)%direction(iptloc) = directemp(jpt) ! store local direction
[2848]	419	ENDIF
[7646]	420	!
	421	END DO
[2848]	422
	423	secs(jsec)%nb_point=iptloc !store number of section's points
	424
	425	!debug
	426	!-----
[8126]	427	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
[2848]	428	WRITE(numout,*)" List of points selected by the proc:"
	429	DO jpt = 1,iptloc
[7646]	430	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
	431	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
[2848]	432	WRITE(numout,*)' # I J : ',iiglo,ijglo
	433	ENDDO
	434	ENDIF
	435
[3294]	436	IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
	437	DO jpt = 1,iptloc
[7646]	438	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
	439	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
[3294]	440	ENDDO
	441	ENDIF
	442
[2848]	443	!remove redundant points between processors
[2908]	444	!------------------------------------------
[8126]	445	lldebug = .FALSE. ; IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) lldebug = .TRUE.
[2848]	446	IF( iptloc .NE. 0 )THEN
	447	CALL removepoints(secs(jsec),'I','top_list',lldebug)
	448	CALL removepoints(secs(jsec),'I','bot_list',lldebug)
	449	CALL removepoints(secs(jsec),'J','top_list',lldebug)
	450	CALL removepoints(secs(jsec),'J','bot_list',lldebug)
	451	ENDIF
[3294]	452	IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
	453	DO jpt = 1,secs(jsec)%nb_point
[7646]	454	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
	455	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
[3294]	456	ENDDO
	457	ENDIF
[2848]	458
	459	!debug
	460	!-----
[8126]	461	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
[2848]	462	WRITE(numout,*)" List of points after removepoints:"
[3294]	463	iptloc = secs(jsec)%nb_point
[2848]	464	DO jpt = 1,iptloc
[7646]	465	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
	466	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
[2848]	467	WRITE(numout,*)' # I J : ',iiglo,ijglo
[4153]	468	CALL FLUSH(numout)
[2848]	469	ENDDO
	470	ENDIF
	471
	472	ELSE ! iptglo = 0
[8126]	473	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )&
[2912]	474	WRITE(numout,*)' No points for this section.'
[2848]	475	ENDIF
	476
	477	ENDDO !end of the loop on jsec
	478
	479	nb_sec = jsec-1 !number of section read in the file
[3294]	480	!
[2848]	481	END SUBROUTINE readsec
	482
[9124]	483
[2848]	484	SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
	485	!!---------------------------------------------------------------------------
	486	!! *** function removepoints
	487	!!
[3680]	488	!! ** Purpose :: Remove points which are common to 2 procs
[2854]	489	!!
[2848]	490	!----------------------------------------------------------------------------
	491	!! * arguments
	492	TYPE(SECTION),INTENT(INOUT) :: sec
	493	CHARACTER(len=1),INTENT(IN) :: cdind ! = 'I'/'J'
	494	CHARACTER(len=8),INTENT(IN) :: cdextr ! = 'top_list'/'bot_list'
	495	LOGICAL,INTENT(IN) :: ld_debug
	496
	497	!! * Local variables
	498	INTEGER :: iextr ,& !extremity of listpoint that we verify
	499	iind ,& !coord of listpoint that we verify
	500	itest ,& !indice value of the side of the domain
	501	!where points could be redundant
[2912]	502	isgn ,& ! isgn= 1 : scan listpoint from start to end
	503	! isgn=-1 : scan listpoint from end to start
[2848]	504	istart,iend !first and last points selected in listpoint
[3294]	505	INTEGER :: jpoint !loop on list points
[9125]	506	INTEGER, DIMENSION(nb_point_max) :: idirec !contains temporary sec%direction
	507	INTEGER, DIMENSION(2,nb_point_max) :: icoord !contains temporary sec%listpoint
[2848]	508	!----------------------------------------------------------------------------
[9125]	509	!
[2848]	510	IF( ld_debug )WRITE(numout,*)' -------------------------'
	511	IF( ld_debug )WRITE(numout,*)' removepoints in listpoint'
	512
	513	!iextr=extremity of list_point that we verify
	514	IF ( cdextr=='bot_list' )THEN ; iextr=1 ; isgn=1
	515	ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
	516	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdextr")
	517	ENDIF
	518
	519	!which coordinate shall we verify ?
	520	IF ( cdind=='I' )THEN ; itest=nlei ; iind=1
	521	ELSE IF ( cdind=='J' )THEN ; itest=nlej ; iind=2
	522	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdind")
	523	ENDIF
	524
	525	IF( ld_debug )THEN
	526	WRITE(numout,*)' case: coord/list extr/domain side'
	527	WRITE(numout,*)' ', cdind,' ',cdextr,' ',itest
	528	WRITE(numout,*)' Actual number of points: ',sec%nb_point
	529	ENDIF
	530
	531	icoord(1,1:nb_point_max) = sec%listPoint%I
	532	icoord(2,1:nb_point_max) = sec%listPoint%J
	533	idirec = sec%direction
	534	sec%listPoint = POINT_SECTION(0,0)
	535	sec%direction = 0
	536
	537	jpoint=iextr+isgn
[3294]	538	DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
	539	IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
	540	ELSE ; EXIT
	541	ENDIF
	542	ENDDO
[2908]	543
[2848]	544	IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
	545	ELSE ; istart=1 ; iend=jpoint+1
	546	ENDIF
[3294]	547
[2848]	548	sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
	549	sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
	550	sec%direction(1:1+iend-istart) = idirec(istart:iend)
	551	sec%nb_point = iend-istart+1
	552
	553	IF( ld_debug )THEN
	554	WRITE(numout,*)' Number of points after removepoints :',sec%nb_point
	555	WRITE(numout,*)' sec%direction after removepoints :',sec%direction(1:sec%nb_point)
	556	ENDIF
[9124]	557	!
	558	END SUBROUTINE removepoints
[2848]	559
[9125]	560
[10965]	561	SUBROUTINE transport(Kmm,sec,ld_debug,jsec)
[2912]	562	!!-------------------------------------------------------------------------------------------
[2848]	563	!! * ROUTINE transport *
	564	!!
[3680]	565	!! Purpose :: Compute the transport for each point in a section
[2913]	566	!!
[3680]	567	!! Method :: Loop over each segment, and each vertical level and add the transport
	568	!! Be aware :
	569	!! One section is a sum of segments
	570	!! One segment is defined by 2 consecutive points in sec%listPoint
	571	!! All points of sec%listPoint are positioned on the F-point of the cell
	572	!!
	573	!! There are two loops:
	574	!! loop on the segment between 2 nodes
	575	!! loop on the level jk !!
	576	!!
	577	!! Output :: Arrays containing the volume,density,heat,salt transports for each i
	578	!! point in a section, summed over each nn_dct.
[2848]	579	!!
[2912]	580	!!-------------------------------------------------------------------------------------------
[10965]	581	INTEGER ,INTENT(IN) :: Kmm ! time level index
[2848]	582	TYPE(SECTION),INTENT(INOUT) :: sec
	583	LOGICAL ,INTENT(IN) :: ld_debug
[3680]	584	INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
[9124]	585	!
	586	INTEGER :: jk, jseg, jclass,jl, isgnu, isgnv ! loop on level/segment/classes/ice categories
	587	REAL(wp):: zumid, zvmid, zumid_ice, zvmid_ice ! U/V ocean & ice velocity on a cell segment
	588	REAL(wp):: zTnorm ! transport of velocity through one cell's sides
	589	REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/potential density/ssh/depth at u/v point
	590	TYPE(POINT_SECTION) :: k
[6140]	591	!!--------------------------------------------------------
	592	!
	593	IF( ld_debug )WRITE(numout,*)' Compute transport'
[2848]	594
[6140]	595	!---------------------------!
	596	! COMPUTE TRANSPORT !
	597	!---------------------------!
	598	IF(sec%nb_point .NE. 0)THEN
[2848]	599
[6140]	600	!----------------------------------------------------------------------------------------------------
	601	!Compute sign for velocities:
	602	!
	603	!convention:
	604	! non horizontal section: direction + is toward left hand of section
	605	! horizontal section: direction + is toward north of section
	606	!
	607	!
	608	! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0
	609	! ---------------- ----------------- --------------- --------------
	610	!
	611	! isgnv=1 direction +
	612	! ______ _____ ______
	613	! \| //\| \| \| direction +
	614	! \| isgnu=1 // \| \|isgnu=1 \|isgnu=1 /\|\
	615	! \|_______ // ______\| \\ \| ---\ \|
	616	! \| \| isgnv=-1 \\ \| \| ---/ direction + ____________
	617	! \| \| __\\\| \|
	618	! \| \| direction + \| isgnv=1
	619	!
	620	!----------------------------------------------------------------------------------------------------
	621	isgnu = 1
	622	IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1
	623	ELSE ; isgnv = 1
	624	ENDIF
	625	IF( sec%slopeSection .GE. 9999. ) isgnv = 1
[2848]	626
[6140]	627	IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
[2848]	628
[6140]	629	!--------------------------------------!
	630	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
	631	!--------------------------------------!
	632	DO jseg=1,MAX(sec%nb_point-1,0)
[2848]	633
[6140]	634	!-------------------------------------------------------------------------------------------
	635	! Select the appropriate coordinate for computing the velocity of the segment
	636	!
	637	! CASE(0) Case (2)
	638	! ------- --------
	639	! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
	640	! F(i,j)----------V(i+1,j)-------F(i+1,j) \|
	641	! \|
	642	! \|
	643	! \|
	644	! Case (3) U(i,j)
	645	! -------- \|
	646	! \|
	647	! listPoint(jseg+1) F(i,j+1) \|
	648	! \| \|
	649	! \| \|
	650	! \| listPoint(jseg+1) F(i,j-1)
	651	! \|
	652	! \|
	653	! U(i,j+1)
	654	! \| Case(1)
	655	! \| ------
	656	! \|
	657	! \| listPoint(jseg+1) listPoint(jseg)
	658	! \| F(i-1,j)-----------V(i,j) -------f(jseg)
	659	! listPoint(jseg) F(i,j)
	660	!
	661	!-------------------------------------------------------------------------------------------
[2919]	662
[6140]	663	SELECT CASE( sec%direction(jseg) )
	664	CASE(0) ; k = sec%listPoint(jseg)
	665	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
	666	CASE(2) ; k = sec%listPoint(jseg)
	667	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
	668	END SELECT
[2864]	669
[6140]	670	!---------------------------\|
	671	! LOOP ON THE LEVEL \|
	672	!---------------------------\|
[8126]	673	DO jk = 1, mbkt(k%I,k%J) !Sum of the transport on the vertical
[6140]	674	! ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
	675	SELECT CASE( sec%direction(jseg) )
	676	CASE(0,1)
[10965]	677	ztn = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) )
	678	zsn = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) )
	679	zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop)
	680	zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0)
	681	zsshn = 0.5( ssh(k%I,k%J,Kmm) + ssh(k%I,k%J+1,Kmm) ) vmask(k%I,k%J,1)
[6140]	682	CASE(2,3)
[10965]	683	ztn = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) )
	684	zsn = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) )
	685	zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop)
	686	zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0)
	687	zsshn = 0.5( ssh(k%I,k%J,Kmm) + ssh(k%I+1,k%J,Kmm) ) umask(k%I,k%J,1)
[6140]	688	END SELECT
	689	!
[10965]	690	zdep= gdept(k%I,k%J,jk,Kmm)
[3680]	691
[6140]	692	SELECT CASE( sec%direction(jseg) ) !compute velocity with the correct direction
	693	CASE(0,1)
	694	zumid=0._wp
[10965]	695	zvmid=isgnvvv(k%I,k%J,jk,Kmm)vmask(k%I,k%J,jk)
[6140]	696	CASE(2,3)
[10965]	697	zumid=isgnuuu(k%I,k%J,jk,Kmm)umask(k%I,k%J,jk)
[6140]	698	zvmid=0._wp
	699	END SELECT
[3680]	700
[6140]	701	!zTnorm=transport through one cell;
	702	!velocity* cell's length * cell's thickness
[10965]	703	zTnorm = zumide2u(k%I,k%J) e3u(k%I,k%J,jk,Kmm) &
	704	& + zvmide1v(k%I,k%J) e3v(k%I,k%J,jk,Kmm)
[2848]	705
[6140]	706	!!gm THIS is WRONG no transport due to ssh in linear free surface case !!!!!
	707	IF( ln_linssh ) THEN !add transport due to free surface
	708	IF( jk==1 ) THEN
	709	zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) &
	710	& + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
	711	ENDIF
	712	ENDIF
	713	!!gm end
[3680]	714	!COMPUTE TRANSPORT
[2848]	715
[3680]	716	transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
	717
[6140]	718	IF( sec%llstrpond ) THEN
[3680]	719	transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * ztn * zrhop * rcp
	720	transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zsn * zrhop * 0.001
[2848]	721	ENDIF
	722
[6140]	723	END DO !end of loop on the level
[2848]	724
[9570]	725	#if defined key_si3
[2848]	726
	727	!ICE CASE
	728	!------------
	729	IF( sec%ll_ice_section )THEN
	730	SELECT CASE (sec%direction(jseg))
	731	CASE(0)
	732	zumid_ice = 0
	733	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
	734	CASE(1)
	735	zumid_ice = 0
	736	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
	737	CASE(2)
	738	zvmid_ice = 0
	739	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
	740	CASE(3)
	741	zvmid_ice = 0
	742	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
	743	END SELECT
	744
	745	zTnorm=zumid_icee2u(k%I,k%J)+zvmid_icee1v(k%I,k%J)
[4153]	746
[9570]	747	#if defined key_si3
[4153]	748	DO jl=1,jpl
[9019]	749	transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* &
	750	a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * &
	751	( h_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) + &
	752	h_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
[4153]	753
	754	transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* &
[9019]	755	a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
[6140]	756	END DO
[4153]	757	#endif
[2848]	758
	759	ENDIF !end of ice case
	760	#endif
	761
[6140]	762	END DO !end of loop on the segment
[2848]	763
[3680]	764	ENDIF !end of sec%nb_point =0 case
[3294]	765	!
[2848]	766	END SUBROUTINE transport
[6140]	767
	768
[10965]	769	SUBROUTINE dia_dct_sum(Kmm,sec,jsec)
[3680]	770	!!-------------------------------------------------------------
	771	!! Purpose: Average the transport over nn_dctwri time steps
	772	!! and sum over the density/salinity/temperature/depth classes
	773	!!
	774	!! Method: Sum over relevant grid cells to obtain values
	775	!! for each class
	776	!! There are several loops:
	777	!! loop on the segment between 2 nodes
	778	!! loop on the level jk
	779	!! loop on the density/temperature/salinity/level classes
	780	!! test on the density/temperature/salinity/level
	781	!!
	782	!! Note: Transport through a given section is equal to the sum of transports
	783	!! computed on each proc.
	784	!! On each proc,transport is equal to the sum of transport computed through
	785	!! segments linking each point of sec%listPoint with the next one.
	786	!!
	787	!!-------------------------------------------------------------
[10965]	788	INTEGER ,INTENT(IN) :: Kmm ! time level index
[3680]	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)
[10965]	850	ztn = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) )
	851	zsn = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) )
	852	zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop)
	853	zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0)
[3680]	854
	855	CASE(2,3)
[10965]	856	ztn = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) )
	857	zsn = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) )
	858	zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop)
	859	zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0)
	860	zsshn = 0.5( ssh(k%I,k%J,Kmm) + ssh(k%I+1,k%J,Kmm) ) umask(k%I,k%J,1)
[3680]	861	END SELECT
	862
[10965]	863	zdep= gdept(k%I,k%J,jk,Kmm)
[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
[10965]	1106	FUNCTION interp(Kmm, 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
[10965]	1167	INTEGER, INTENT(IN) :: Kmm ! time level index
[2908]	1168	INTEGER, INTENT(IN) :: ki, kj, kk ! coordinate of point
	1169	CHARACTER(len=1), INTENT(IN) :: cd_point ! type of point (U, V)
	1170	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab ! variable to compute at (ki, kj, kk )
	1171	REAL(wp) :: interp ! interpolated variable
[2848]	1172
	1173	!*local declations
[2908]	1174	INTEGER :: ii1, ij1, ii2, ij2 ! local integer
[6140]	1175	REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu ! local real
[2908]	1176	REAL(wp):: zet1, zet2 ! weight for interpolation
	1177	REAL(wp):: zdep1,zdep2 ! differences of depth
[4613]	1178	REAL(wp):: zmsk ! mask value
[2848]	1179	!!----------------------------------------------------------------------
	1180
	1181	IF( cd_point=='U' )THEN
	1182	ii1 = ki ; ij1 = kj
	1183	ii2 = ki+1 ; ij2 = kj
[2908]	1184
	1185	zet1=e1t(ii1,ij1)
	1186	zet2=e1t(ii2,ij2)
[4613]	1187	zmsk=umask(ii1,ij1,kk)
	1188
[2908]	1189
[2848]	1190	ELSE ! cd_point=='V'
	1191	ii1 = ki ; ij1 = kj
	1192	ii2 = ki ; ij2 = kj+1
[2908]	1193
	1194	zet1=e2t(ii1,ij1)
	1195	zet2=e2t(ii2,ij2)
[4613]	1196	zmsk=vmask(ii1,ij1,kk)
[2908]	1197
[2848]	1198	ENDIF
	1199
[2908]	1200	IF( ln_sco )THEN ! s-coordinate case
	1201
[10965]	1202	zdepu = ( gdept(ii1,ij1,kk,Kmm) + gdept(ii2,ij2,kk,Kmm) ) * 0.5_wp
	1203	zdep1 = gdept(ii1,ij1,kk,Kmm) - zdepu
	1204	zdep2 = gdept(ii2,ij2,kk,Kmm) - zdepu
[2908]	1205
[3680]	1206	! weights
[2908]	1207	zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
	1208	zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
	1209
	1210	! result
[4613]	1211	interp = zmsk * ( zwgt2 * ptab(ii1,ij1,kk) + zwgt1 * ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )
[2908]	1212
	1213
	1214	ELSE ! full step or partial step case
	1215
[10965]	1216	ze3t = e3t(ii2,ij2,kk,Kmm) - e3t(ii1,ij1,kk,Kmm)
	1217	zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) ) / e3w(ii2,ij2,kk,Kmm)
	1218	zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) ) / e3w(ii1,ij1,kk,Kmm)
[2848]	1219
[2908]	1220	IF(kk .NE. 1)THEN
[2848]	1221
[2908]	1222	IF( ze3t >= 0. )THEN
[3680]	1223	! zbis
[2908]	1224	zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )
	1225	! result
[4613]	1226	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
[2908]	1227	ELSE
[3680]	1228	! zbis
[2908]	1229	zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
	1230	! result
[4613]	1231	interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
[2908]	1232	ENDIF
	1233
[2848]	1234	ELSE
[4613]	1235	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
[2908]	1236	ENDIF
[2848]	1237
	1238	ENDIF
[6140]	1239	!
	1240	END FUNCTION interp
[2848]	1241
	1242	#else
	1243	!!----------------------------------------------------------------------
[11822]	1244	!! Dummy module
[2848]	1245	!!----------------------------------------------------------------------
[11822]	1246	LOGICAL, PUBLIC :: ln_diadct = .FALSE.
[2848]	1247	CONTAINS
[11822]	1248	SUBROUTINE dia_dct_init
[9927]	1249	IMPLICIT NONE
[2864]	1250	END SUBROUTINE dia_dct_init
	1251
[10965]	1252	SUBROUTINE dia_dct( kt, Kmm ) ! Dummy routine
[9927]	1253	IMPLICIT NONE
[3680]	1254	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[10965]	1255	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
[2848]	1256	WRITE(,) 'dia_dct: You should not have seen this print! error?', kt
	1257	END SUBROUTINE dia_dct
[11822]	1258	!
[2848]	1259	#endif
	1260
[6140]	1261	!!======================================================================
[2848]	1262	END MODULE diadct

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

Download in other formats: