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_r3342_MERCATOR7_SST/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2012/dev_r3342_MERCATOR7_SST/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 3572

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

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

Download in other formats: