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 branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 3988

Last change on this file since 3988 was 3988, checked in by cbricaud, 11 years ago
bugfix diadct for ice transport ( LIM3 case ) ; see ticket 1080
Property svn:executable set to ``*
File size: 54.9 KB

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

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

Download in other formats: