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

Context Navigation

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 4147

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

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

Download in other formats: