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_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diadct.F90 @ 12178

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

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

Download in other formats: