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

Context Navigation

source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 5726

Last change on this file since 5726 was 5726, checked in by jpalmier, 9 years ago
JPALM -- 10-09-2015 -- add MEDUSA in the branch ; adapted TOP_SRC to MEDUSA ; remove some svn keywords in the branch
Property svn:executable set to ``*
File size: 59.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	!!----------------------------------------------------------------------
[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
[4153]	44	USE ice
[3632]	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 )
[4624]	152	IF(lwm) WRITE ( numond, namdct )
[4147]	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
[5505]	177	IF( lwm ) THEN
[2927]	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
[5505]	284	IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
[2848]	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
[4153]	492	CALL FLUSH(numout)
[2848]	493	ENDDO
	494	ENDIF
	495
	496	ELSE ! iptglo = 0
	497	IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )&
[2912]	498	WRITE(numout,*)' No points for this section.'
[2848]	499	ENDIF
	500
	501	ENDDO !end of the loop on jsec
	502
	503	nb_sec = jsec-1 !number of section read in the file
	504
[3294]	505	CALL wrk_dealloc( nb_point_max, directemp )
	506	!
[2848]	507	END SUBROUTINE readsec
	508
	509	SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
	510	!!---------------------------------------------------------------------------
	511	!! *** function removepoints
	512	!!
[3680]	513	!! ** Purpose :: Remove points which are common to 2 procs
[2854]	514	!!
[2848]	515	!----------------------------------------------------------------------------
	516	!! * arguments
	517	TYPE(SECTION),INTENT(INOUT) :: sec
	518	CHARACTER(len=1),INTENT(IN) :: cdind ! = 'I'/'J'
	519	CHARACTER(len=8),INTENT(IN) :: cdextr ! = 'top_list'/'bot_list'
	520	LOGICAL,INTENT(IN) :: ld_debug
	521
	522	!! * Local variables
	523	INTEGER :: iextr ,& !extremity of listpoint that we verify
	524	iind ,& !coord of listpoint that we verify
	525	itest ,& !indice value of the side of the domain
	526	!where points could be redundant
[2912]	527	isgn ,& ! isgn= 1 : scan listpoint from start to end
	528	! isgn=-1 : scan listpoint from end to start
[2848]	529	istart,iend !first and last points selected in listpoint
[3294]	530	INTEGER :: jpoint !loop on list points
	531	INTEGER, POINTER, DIMENSION(:) :: idirec !contains temporary sec%direction
	532	INTEGER, POINTER, DIMENSION(:,:) :: icoord !contains temporary sec%listpoint
[2848]	533	!----------------------------------------------------------------------------
[3294]	534	CALL wrk_alloc( nb_point_max, idirec )
	535	CALL wrk_alloc( 2, nb_point_max, icoord )
	536
[2848]	537	IF( ld_debug )WRITE(numout,*)' -------------------------'
	538	IF( ld_debug )WRITE(numout,*)' removepoints in listpoint'
	539
	540	!iextr=extremity of list_point that we verify
	541	IF ( cdextr=='bot_list' )THEN ; iextr=1 ; isgn=1
	542	ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
	543	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdextr")
	544	ENDIF
	545
	546	!which coordinate shall we verify ?
	547	IF ( cdind=='I' )THEN ; itest=nlei ; iind=1
	548	ELSE IF ( cdind=='J' )THEN ; itest=nlej ; iind=2
	549	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdind")
	550	ENDIF
	551
	552	IF( ld_debug )THEN
	553	WRITE(numout,*)' case: coord/list extr/domain side'
	554	WRITE(numout,*)' ', cdind,' ',cdextr,' ',itest
	555	WRITE(numout,*)' Actual number of points: ',sec%nb_point
	556	ENDIF
	557
	558	icoord(1,1:nb_point_max) = sec%listPoint%I
	559	icoord(2,1:nb_point_max) = sec%listPoint%J
	560	idirec = sec%direction
	561	sec%listPoint = POINT_SECTION(0,0)
	562	sec%direction = 0
	563
	564	jpoint=iextr+isgn
[3294]	565	DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
	566	IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
	567	ELSE ; EXIT
	568	ENDIF
	569	ENDDO
[2908]	570
[2848]	571	IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
	572	ELSE ; istart=1 ; iend=jpoint+1
	573	ENDIF
[3294]	574
[2848]	575	sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
	576	sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
	577	sec%direction(1:1+iend-istart) = idirec(istart:iend)
	578	sec%nb_point = iend-istart+1
	579
	580	IF( ld_debug )THEN
	581	WRITE(numout,*)' Number of points after removepoints :',sec%nb_point
	582	WRITE(numout,*)' sec%direction after removepoints :',sec%direction(1:sec%nb_point)
	583	ENDIF
	584
[3294]	585	CALL wrk_dealloc( nb_point_max, idirec )
	586	CALL wrk_dealloc( 2, nb_point_max, icoord )
[2848]	587	END SUBROUTINE removepoints
	588
[3680]	589	SUBROUTINE transport(sec,ld_debug,jsec)
[2912]	590	!!-------------------------------------------------------------------------------------------
[2848]	591	!! * ROUTINE transport *
	592	!!
[3680]	593	!! Purpose :: Compute the transport for each point in a section
[2913]	594	!!
[3680]	595	!! Method :: Loop over each segment, and each vertical level and add the transport
	596	!! Be aware :
	597	!! One section is a sum of segments
	598	!! One segment is defined by 2 consecutive points in sec%listPoint
	599	!! All points of sec%listPoint are positioned on the F-point of the cell
	600	!!
	601	!! There are two loops:
	602	!! loop on the segment between 2 nodes
	603	!! loop on the level jk !!
	604	!!
	605	!! Output :: Arrays containing the volume,density,heat,salt transports for each i
	606	!! point in a section, summed over each nn_dct.
[2848]	607	!!
[2912]	608	!!-------------------------------------------------------------------------------------------
[2848]	609	!! * Arguments
	610	TYPE(SECTION),INTENT(INOUT) :: sec
	611	LOGICAL ,INTENT(IN) :: ld_debug
[3680]	612	INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
[2848]	613
	614	!! * Local variables
[4153]	615	INTEGER :: jk, jseg, jclass,jl, &!loop on level/segment/classes/ice categories
[3680]	616	isgnu, isgnv !
	617	REAL(wp) :: zumid, zvmid, &!U/V velocity on a cell segment
	618	zumid_ice, zvmid_ice, &!U/V ice velocity
	619	zTnorm !transport of velocity through one cell's sides
	620	REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep !temperature/salinity/potential density/ssh/depth at u/v point
[2848]	621
	622	TYPE(POINT_SECTION) :: k
	623	!!--------------------------------------------------------
	624
	625	IF( ld_debug )WRITE(numout,*)' Compute transport'
	626
	627	!---------------------------!
	628	! COMPUTE TRANSPORT !
	629	!---------------------------!
	630	IF(sec%nb_point .NE. 0)THEN
	631
	632	!----------------------------------------------------------------------------------------------------
	633	!Compute sign for velocities:
	634	!
	635	!convention:
[2912]	636	! non horizontal section: direction + is toward left hand of section
	637	! horizontal section: direction + is toward north of section
[2848]	638	!
	639	!
	640	! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0
	641	! ---------------- ----------------- --------------- --------------
	642	!
[2912]	643	! isgnv=1 direction +
	644	! ______ _____ ______
	645	! \| //\| \| \| direction +
	646	! \| isgnu=1 // \| \|isgnu=1 \|isgnu=1 /\|\
	647	! \|_______ // ______\| \\ \| ---\ \|
	648	! \| \| isgnv=-1 \\ \| \| ---/ direction + ____________
	649	! \| \| __\\\| \|
	650	! \| \| direction + \| isgnv=1
	651	!
[2848]	652	!----------------------------------------------------------------------------------------------------
	653	isgnu = 1
	654	IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1
	655	ELSE ; isgnv = 1
	656	ENDIF
[3632]	657	IF( sec%slopeSection .GE. 9999. ) isgnv = 1
[2848]	658
[3632]	659	IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
[2848]	660
	661	!--------------------------------------!
	662	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
	663	!--------------------------------------!
	664	DO jseg=1,MAX(sec%nb_point-1,0)
	665
[2919]	666	!-------------------------------------------------------------------------------------------
[2927]	667	! Select the appropriate coordinate for computing the velocity of the segment
[2848]	668	!
[2919]	669	! CASE(0) Case (2)
	670	! ------- --------
	671	! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
	672	! F(i,j)----------V(i+1,j)-------F(i+1,j) \|
	673	! \|
	674	! \|
	675	! \|
	676	! Case (3) U(i,j)
	677	! -------- \|
	678	! \|
	679	! listPoint(jseg+1) F(i,j+1) \|
	680	! \| \|
	681	! \| \|
	682	! \| listPoint(jseg+1) F(i,j-1)
	683	! \|
	684	! \|
	685	! U(i,j+1)
	686	! \| Case(1)
	687	! \| ------
	688	! \|
	689	! \| listPoint(jseg+1) listPoint(jseg)
	690	! \| F(i-1,j)-----------V(i,j) -------f(jseg)
	691	! listPoint(jseg) F(i,j)
	692	!
	693	!-------------------------------------------------------------------------------------------
	694
[2848]	695	SELECT CASE( sec%direction(jseg) )
	696	CASE(0) ; k = sec%listPoint(jseg)
	697	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
	698	CASE(2) ; k = sec%listPoint(jseg)
	699	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
	700	END SELECT
[2864]	701
[3680]	702	!---------------------------\|
	703	! LOOP ON THE LEVEL \|
	704	!---------------------------\|
	705	!Sum of the transport on the vertical
	706	DO jk=1,mbathy(k%I,k%J)
[2848]	707
[3680]	708	! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
	709	SELECT CASE( sec%direction(jseg) )
	710	CASE(0,1)
	711	ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
	712	zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
	713	zrhop = interp(k%I,k%J,jk,'V',rhop)
	714	zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
	715	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) vmask(k%I,k%J,1)
	716	CASE(2,3)
	717	ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
	718	zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
	719	zrhop = interp(k%I,k%J,jk,'U',rhop)
	720	zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
	721	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) umask(k%I,k%J,1)
	722	END SELECT
[2848]	723
[4613]	724	zfsdep= fsdept(k%I,k%J,jk)
[3680]	725
	726	!compute velocity with the correct direction
	727	SELECT CASE( sec%direction(jseg) )
	728	CASE(0,1)
	729	zumid=0.
	730	zvmid=isgnvvn(k%I,k%J,jk)vmask(k%I,k%J,jk)
	731	CASE(2,3)
	732	zumid=isgnuun(k%I,k%J,jk)umask(k%I,k%J,jk)
	733	zvmid=0.
	734	END SELECT
	735
	736	!zTnorm=transport through one cell;
	737	!velocity* cell's length * cell's thickness
	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)
[2848]	740
	741	#if ! defined key_vvl
[3680]	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
[2848]	747	#endif
[3680]	748	!COMPUTE TRANSPORT
[2848]	749
[3680]	750	transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
	751
	752	IF ( sec%llstrpond ) THEN
	753	transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * ztn * zrhop * rcp
	754	transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zsn * zrhop * 0.001
[2848]	755	ENDIF
	756
[3680]	757	ENDDO !end of loop on the level
[2848]	758
[3294]	759	#if defined key_lim2 \|\| defined key_lim3
[2848]	760
	761	!ICE CASE
	762	!------------
	763	IF( sec%ll_ice_section )THEN
	764	SELECT CASE (sec%direction(jseg))
	765	CASE(0)
	766	zumid_ice = 0
	767	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
	768	CASE(1)
	769	zumid_ice = 0
	770	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
	771	CASE(2)
	772	zvmid_ice = 0
	773	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
	774	CASE(3)
	775	zvmid_ice = 0
	776	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
	777	END SELECT
	778
	779	zTnorm=zumid_icee2u(k%I,k%J)+zvmid_icee1v(k%I,k%J)
[4153]	780
	781	#if defined key_lim2
[3680]	782	transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* &
	783	(1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
	784	*(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + &
	785	hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
	786	transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* &
	787	(1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
[4153]	788	#endif
	789	#if defined key_lim3
	790	DO jl=1,jpl
	791	transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* &
	792	a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * &
	793	( ht_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) + &
	794	ht_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
	795
	796	transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* &
	797	a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
	798	ENDDO
	799	#endif
[2848]	800
	801	ENDIF !end of ice case
	802	#endif
	803
	804	ENDDO !end of loop on the segment
	805
[3680]	806	ENDIF !end of sec%nb_point =0 case
[3294]	807	!
[2848]	808	END SUBROUTINE transport
	809
[3680]	810	SUBROUTINE dia_dct_sum(sec,jsec)
	811	!!-------------------------------------------------------------
	812	!! Purpose: Average the transport over nn_dctwri time steps
	813	!! and sum over the density/salinity/temperature/depth classes
	814	!!
	815	!! Method: Sum over relevant grid cells to obtain values
	816	!! for each class
	817	!! There are several loops:
	818	!! loop on the segment between 2 nodes
	819	!! loop on the level jk
	820	!! loop on the density/temperature/salinity/level classes
	821	!! test on the density/temperature/salinity/level
	822	!!
	823	!! Note: Transport through a given section is equal to the sum of transports
	824	!! computed on each proc.
	825	!! On each proc,transport is equal to the sum of transport computed through
	826	!! segments linking each point of sec%listPoint with the next one.
	827	!!
	828	!!-------------------------------------------------------------
	829	!! * arguments
	830	TYPE(SECTION),INTENT(INOUT) :: sec
	831	INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
	832
	833	TYPE(POINT_SECTION) :: k
	834	INTEGER :: jk,jseg,jclass ! dummy variables for looping on level/segment/classes
	835	REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
	836	!!-------------------------------------------------------------
	837
	838	!! Sum the relevant segments to obtain values for each class
	839	IF(sec%nb_point .NE. 0)THEN
	840
	841	!--------------------------------------!
	842	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
	843	!--------------------------------------!
	844	DO jseg=1,MAX(sec%nb_point-1,0)
	845
	846	!-------------------------------------------------------------------------------------------
	847	! Select the appropriate coordinate for computing the velocity of the segment
	848	!
	849	! CASE(0) Case (2)
	850	! ------- --------
	851	! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
	852	! F(i,j)----------V(i+1,j)-------F(i+1,j) \|
	853	! \|
	854	! \|
	855	! \|
	856	! Case (3) U(i,j)
	857	! -------- \|
	858	! \|
	859	! listPoint(jseg+1) F(i,j+1) \|
	860	! \| \|
	861	! \| \|
	862	! \| listPoint(jseg+1) F(i,j-1)
	863	! \|
	864	! \|
	865	! U(i,j+1)
	866	! \| Case(1)
	867	! \| ------
	868	! \|
	869	! \| listPoint(jseg+1) listPoint(jseg)
	870	! \| F(i-1,j)-----------V(i,j) -------f(jseg)
	871	! listPoint(jseg) F(i,j)
	872	!
	873	!-------------------------------------------------------------------------------------------
	874
	875	SELECT CASE( sec%direction(jseg) )
	876	CASE(0) ; k = sec%listPoint(jseg)
	877	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
	878	CASE(2) ; k = sec%listPoint(jseg)
	879	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
	880	END SELECT
	881
	882	!---------------------------\|
	883	! LOOP ON THE LEVEL \|
	884	!---------------------------\|
	885	!Sum of the transport on the vertical
	886	DO jk=1,mbathy(k%I,k%J)
	887
	888	! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
	889	SELECT CASE( sec%direction(jseg) )
	890	CASE(0,1)
	891	ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
	892	zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
	893	zrhop = interp(k%I,k%J,jk,'V',rhop)
	894	zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
	895
	896	CASE(2,3)
	897	ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
	898	zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
	899	zrhop = interp(k%I,k%J,jk,'U',rhop)
	900	zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
	901	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) umask(k%I,k%J,1)
	902	END SELECT
	903
[4613]	904	zfsdep= fsdept(k%I,k%J,jk)
[3680]	905
	906	!-------------------------------
	907	! LOOP ON THE DENSITY CLASSES \|
	908	!-------------------------------
	909	!The computation is made for each density/temperature/salinity/depth class
	910	DO jclass=1,MAX(1,sec%nb_class-1)
	911
	912	!----------------------------------------------!
	913	!TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
	914	!----------------------------------------------!
	915
	916	IF ( ( &
	917	((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. &
	918	( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. &
	919	( sec%zsigp(jclass) .EQ. 99.)) .AND. &
	920
	921	((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. &
	922	( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. &
	923	( sec%zsigi(jclass) .EQ. 99.)) .AND. &
	924
	925	((( zsn .GT. sec%zsal(jclass)) .AND. &
	926	( zsn .LE. sec%zsal(jclass+1))) .OR. &
	927	( sec%zsal(jclass) .EQ. 99.)) .AND. &
	928
	929	((( ztn .GE. sec%ztem(jclass)) .AND. &
	930	( ztn .LE. sec%ztem(jclass+1))) .OR. &
	931	( sec%ztem(jclass) .EQ.99.)) .AND. &
	932
	933	((( zfsdep .GE. sec%zlay(jclass)) .AND. &
	934	( zfsdep .LE. sec%zlay(jclass+1))) .OR. &
	935	( sec%zlay(jclass) .EQ. 99. )) &
	936	)) THEN
	937
	938	!SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
	939	!----------------------------------------------------------------------------
	940	IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN
	941	sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6
	942	ELSE
	943	sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6
	944	ENDIF
	945	IF( sec%llstrpond )THEN
	946
	947	IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN
	948	sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)
	949	ELSE
	950	sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)
	951	ENDIF
	952
	953	IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN
	954	sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)
	955	ELSE
	956	sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)
	957	ENDIF
	958
	959	ELSE
	960	sec%transport( 3,jclass) = 0._wp
	961	sec%transport( 4,jclass) = 0._wp
	962	sec%transport( 5,jclass) = 0._wp
	963	sec%transport( 6,jclass) = 0._wp
	964	ENDIF
	965
	966	ENDIF ! end of test if point is in class
	967
	968	ENDDO ! end of loop on the classes
	969
	970	ENDDO ! loop over jk
	971
	972	#if defined key_lim2 \|\| defined key_lim3
	973
	974	!ICE CASE
	975	IF( sec%ll_ice_section )THEN
	976
	977	IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
	978	sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6
	979	ELSE
	980	sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6
	981	ENDIF
	982
	983	IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
	984	sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6
	985	ELSE
	986	sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6
	987	ENDIF
	988
	989	ENDIF !end of ice case
	990	#endif
	991
	992	ENDDO !end of loop on the segment
	993
	994	ELSE !if sec%nb_point =0
	995	sec%transport(1:2,:)=0.
	996	IF (sec%llstrpond) sec%transport(3:6,:)=0.
	997	IF (sec%ll_ice_section) sec%transport(7:10,:)=0.
	998	ENDIF !end of sec%nb_point =0 case
	999
	1000	END SUBROUTINE dia_dct_sum
	1001
[2941]	1002	SUBROUTINE dia_dct_wri(kt,ksec,sec)
[2848]	1003	!!-------------------------------------------------------------
	1004	!! Write transport output in numdct
	1005	!!
[2854]	1006	!! Purpose: Write transports in ascii files
	1007	!!
	1008	!! Method:
	1009	!! 1. Write volume transports in "volume_transport"
[2908]	1010	!! Unit: Sv : area * Velocity / 1.e6
[2854]	1011	!!
	1012	!! 2. Write heat transports in "heat_transport"
[3680]	1013	!! Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
[2854]	1014	!!
	1015	!! 3. Write salt transports in "salt_transport"
[3680]	1016	!! Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9
[2854]	1017	!!
[2848]	1018	!!-------------------------------------------------------------
	1019	!!arguments
[3294]	1020	INTEGER, INTENT(IN) :: kt ! time-step
	1021	TYPE(SECTION), INTENT(INOUT) :: sec ! section to write
	1022	INTEGER ,INTENT(IN) :: ksec ! section number
[2848]	1023
	1024	!!local declarations
[3680]	1025	INTEGER :: jclass ! Dummy loop
[3294]	1026	CHARACTER(len=2) :: classe ! Classname
	1027	REAL(wp) :: zbnd1,zbnd2 ! Class bounds
	1028	REAL(wp) :: zslope ! section's slope coeff
	1029	!
[3680]	1030	REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
[2848]	1031	!!-------------------------------------------------------------
[3680]	1032	CALL wrk_alloc(nb_type_class , zsumclasses )
[3294]	1033
[3680]	1034	zsumclasses(:)=0._wp
[2908]	1035	zslope = sec%slopeSection
	1036
	1037
[3680]	1038	DO jclass=1,MAX(1,sec%nb_class-1)
[2848]	1039
	1040	classe = 'N '
	1041	zbnd1 = 0._wp
	1042	zbnd2 = 0._wp
[3680]	1043	zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
[2848]	1044
	1045
	1046	!insitu density classes transports
[3680]	1047	IF( ( sec%zsigi(jclass) .NE. 99._wp ) .AND. &
	1048	( sec%zsigi(jclass+1) .NE. 99._wp ) )THEN
[2848]	1049	classe = 'DI '
[3680]	1050	zbnd1 = sec%zsigi(jclass)
	1051	zbnd2 = sec%zsigi(jclass+1)
[2848]	1052	ENDIF
	1053	!potential density classes transports
[3680]	1054	IF( ( sec%zsigp(jclass) .NE. 99._wp ) .AND. &
	1055	( sec%zsigp(jclass+1) .NE. 99._wp ) )THEN
[2848]	1056	classe = 'DP '
[3680]	1057	zbnd1 = sec%zsigp(jclass)
	1058	zbnd2 = sec%zsigp(jclass+1)
[2848]	1059	ENDIF
	1060	!depth classes transports
[3680]	1061	IF( ( sec%zlay(jclass) .NE. 99._wp ) .AND. &
	1062	( sec%zlay(jclass+1) .NE. 99._wp ) )THEN
[2848]	1063	classe = 'Z '
[3680]	1064	zbnd1 = sec%zlay(jclass)
	1065	zbnd2 = sec%zlay(jclass+1)
[2848]	1066	ENDIF
	1067	!salinity classes transports
[3680]	1068	IF( ( sec%zsal(jclass) .NE. 99._wp ) .AND. &
	1069	( sec%zsal(jclass+1) .NE. 99._wp ) )THEN
[2848]	1070	classe = 'S '
[3680]	1071	zbnd1 = sec%zsal(jclass)
	1072	zbnd2 = sec%zsal(jclass+1)
[2848]	1073	ENDIF
	1074	!temperature classes transports
[3680]	1075	IF( ( sec%ztem(jclass) .NE. 99._wp ) .AND. &
	1076	( sec%ztem(jclass+1) .NE. 99._wp ) ) THEN
[2848]	1077	classe = 'T '
[3680]	1078	zbnd1 = sec%ztem(jclass)
	1079	zbnd2 = sec%ztem(jclass+1)
[2848]	1080	ENDIF
	1081
	1082	!write volume transport per class
[2941]	1083	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
[3680]	1084	jclass,classe,zbnd1,zbnd2,&
	1085	sec%transport(1,jclass),sec%transport(2,jclass), &
	1086	sec%transport(1,jclass)+sec%transport(2,jclass)
[2848]	1087
	1088	IF( sec%llstrpond )THEN
	1089
[2854]	1090	!write heat transport per class:
[2941]	1091	WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
[3680]	1092	jclass,classe,zbnd1,zbnd2,&
	1093	sec%transport(3,jclass)1.e-15,sec%transport(4,jclass)1.e-15, &
	1094	( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15
[2848]	1095	!write salt transport per class
[2941]	1096	WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
[3680]	1097	jclass,classe,zbnd1,zbnd2,&
	1098	sec%transport(5,jclass)1.e-9,sec%transport(6,jclass)1.e-9,&
	1099	(sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9
[2848]	1100	ENDIF
	1101
	1102	ENDDO
	1103
	1104	zbnd1 = 0._wp
	1105	zbnd2 = 0._wp
[3680]	1106	jclass=0
[2848]	1107
	1108	!write total volume transport
[2941]	1109	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
[3680]	1110	jclass,"total",zbnd1,zbnd2,&
	1111	zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
[2848]	1112
	1113	IF( sec%llstrpond )THEN
	1114
	1115	!write total heat transport
[2941]	1116	WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
[3680]	1117	jclass,"total",zbnd1,zbnd2,&
	1118	zsumclasses(3)1.e-15,zsumclasses(4)1.e-15,&
	1119	(zsumclasses(3)+zsumclasses(4) )*1.e-15
[2848]	1120	!write total salt transport
[2941]	1121	WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
[3680]	1122	jclass,"total",zbnd1,zbnd2,&
	1123	zsumclasses(5)1.e-9,zsumclasses(6)1.e-9,&
	1124	(zsumclasses(5)+zsumclasses(6))*1.e-9
[2848]	1125	ENDIF
	1126
	1127
	1128	IF ( sec%ll_ice_section) THEN
	1129	!write total ice volume transport
[2941]	1130	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
[3680]	1131	jclass,"ice_vol",zbnd1,zbnd2,&
	1132	sec%transport(7,1),sec%transport(8,1),&
	1133	sec%transport(7,1)+sec%transport(8,1)
[2848]	1134	!write total ice surface transport
[2941]	1135	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
[3680]	1136	jclass,"ice_surf",zbnd1,zbnd2,&
	1137	sec%transport(9,1),sec%transport(10,1), &
	1138	sec%transport(9,1)+sec%transport(10,1)
[2848]	1139	ENDIF
	1140
	1141	118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
	1142	119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
	1143
[3680]	1144	CALL wrk_dealloc(nb_type_class , zsumclasses )
[2848]	1145	END SUBROUTINE dia_dct_wri
	1146
	1147	FUNCTION interp(ki, kj, kk, cd_point, ptab)
	1148	!!----------------------------------------------------------------------
	1149	!!
[3680]	1150	!! Purpose: compute temperature/salinity/density at U-point or V-point
[2908]	1151	!! --------
[2848]	1152	!!
[2908]	1153	!! Method:
	1154	!! ------
[2848]	1155	!!
[2908]	1156	!! ====> full step and partial step
	1157	!!
	1158	!!
[3680]	1159	!! \| I \| I+1 \| Z=temperature/salinity/density at U-poinT
[2908]	1160	!! \| \| \|
[3680]	1161	!! ---------------------------------------- 1. Veritcal interpolation: compute zbis
[2908]	1162	!! \| \| \| interpolation between ptab(I,J,K) and ptab(I,J,K+1)
	1163	!! \| \| \| zbis =
	1164	!! \| \| \| [ e3w(I+1,J,K)ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) ptab(I,J,K-1) ]
	1165	!! \| \| \| /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
	1166	!! \| \| \|
	1167	!! \| \| \| 2. Horizontal interpolation: compute value at U/V point
	1168	!!K-1 \| ptab(I,J,K-1) \| \| interpolation between zbis and ptab(I+1,J,K)
	1169	!! \| . \| \|
[2909]	1170	!! \| . \| \| interp = ( 0.5zet2zbis + 0.5zet1ptab(I+1,J,K) )/(0.5zet2+0.5zet1)
[2908]	1171	!! \| . \| \|
	1172	!! ------------------------------------------
	1173	!! \| . \| \|
	1174	!! \| . \| \|
	1175	!! \| . \| \|
[2909]	1176	!!K \| zbis.......U...ptab(I+1,J,K) \|
[2908]	1177	!! \| . \| \|
	1178	!! \| ptab(I,J,K) \| \|
	1179	!! \| \|------------------\|
	1180	!! \| \| partials \|
	1181	!! \| \| steps \|
	1182	!! -------------------------------------------
	1183	!! <----zet1------><----zet2--------->
	1184	!!
	1185	!!
	1186	!! ====> s-coordinate
	1187	!!
[2909]	1188	!! \| \| \| 1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
	1189	!! \| \| \| Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
[2908]	1190	!! \| \| ptab(I+1,J,K) \|
	1191	!! \| \| T2 \| 2. Interpolation between T1 and T2 values at U point
	1192	!! \| \| ^ \|
	1193	!! \| \| \| zdep2 \|
	1194	!! \| \| \| \|
	1195	!! \| ^ U v \|
	1196	!! \| \| \| \|
	1197	!! \| \| zdep1 \| \|
	1198	!! \| v \| \|
	1199	!! \| T1 \| \|
	1200	!! \| ptab(I,J,K) \| \|
	1201	!! \| \| \|
	1202	!! \| \| \|
	1203	!!
	1204	!! <----zet1--------><----zet2--------->
	1205	!!
[2848]	1206	!!----------------------------------------------------------------------
	1207	!*arguments
[2908]	1208	INTEGER, INTENT(IN) :: ki, kj, kk ! coordinate of point
	1209	CHARACTER(len=1), INTENT(IN) :: cd_point ! type of point (U, V)
	1210	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab ! variable to compute at (ki, kj, kk )
	1211	REAL(wp) :: interp ! interpolated variable
[2848]	1212
	1213	!*local declations
[2908]	1214	INTEGER :: ii1, ij1, ii2, ij2 ! local integer
	1215	REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu ! local real
	1216	REAL(wp):: zet1, zet2 ! weight for interpolation
	1217	REAL(wp):: zdep1,zdep2 ! differences of depth
[4613]	1218	REAL(wp):: zmsk ! mask value
[2848]	1219	!!----------------------------------------------------------------------
	1220
	1221	IF( cd_point=='U' )THEN
	1222	ii1 = ki ; ij1 = kj
	1223	ii2 = ki+1 ; ij2 = kj
[2908]	1224
	1225	zet1=e1t(ii1,ij1)
	1226	zet2=e1t(ii2,ij2)
[4613]	1227	zmsk=umask(ii1,ij1,kk)
	1228
[2908]	1229
[2848]	1230	ELSE ! cd_point=='V'
	1231	ii1 = ki ; ij1 = kj
	1232	ii2 = ki ; ij2 = kj+1
[2908]	1233
	1234	zet1=e2t(ii1,ij1)
	1235	zet2=e2t(ii2,ij2)
[4613]	1236	zmsk=vmask(ii1,ij1,kk)
[2908]	1237
[2848]	1238	ENDIF
	1239
[2908]	1240	IF( ln_sco )THEN ! s-coordinate case
	1241
	1242	zdepu = ( fsdept(ii1,ij1,kk) + fsdept(ii2,ij2,kk) ) /2
	1243	zdep1 = fsdept(ii1,ij1,kk) - zdepu
	1244	zdep2 = fsdept(ii2,ij2,kk) - zdepu
	1245
[3680]	1246	! weights
[2908]	1247	zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
	1248	zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
	1249
	1250	! result
[4613]	1251	interp = zmsk * ( zwgt2 * ptab(ii1,ij1,kk) + zwgt1 * ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )
[2908]	1252
	1253
	1254	ELSE ! full step or partial step case
	1255
[2848]	1256	#if defined key_vvl
	1257
[2927]	1258	ze3t = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk)
	1259	zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk)
	1260	zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk)
[2848]	1261
	1262	#else
	1263
[2927]	1264	ze3t = fse3t(ii2,ij2,kk) - fse3t(ii1,ij1,kk)
[2908]	1265	zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
	1266	zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
[2848]	1267
	1268	#endif
	1269
[2908]	1270	IF(kk .NE. 1)THEN
[2848]	1271
[2908]	1272	IF( ze3t >= 0. )THEN
[3680]	1273	! zbis
[2908]	1274	zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )
	1275	! result
[4613]	1276	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
[2908]	1277	ELSE
[3680]	1278	! zbis
[2908]	1279	zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
	1280	! result
[4613]	1281	interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
[2908]	1282	ENDIF
	1283
[2848]	1284	ELSE
[4613]	1285	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
[2908]	1286	ENDIF
[2848]	1287
	1288	ENDIF
	1289
[2908]	1290
[2848]	1291	END FUNCTION interp
	1292
	1293	#else
	1294	!!----------------------------------------------------------------------
	1295	!! Default option : Dummy module
	1296	!!----------------------------------------------------------------------
	1297	LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .FALSE. !: diamht flag
[2864]	1298	PUBLIC
[2848]	1299	CONTAINS
[2864]	1300
	1301	SUBROUTINE dia_dct_init ! Dummy routine
	1302	WRITE(,) 'dia_dct_init: You should not have seen this print! error?', kt
	1303	END SUBROUTINE dia_dct_init
	1304
[3680]	1305	SUBROUTINE dia_dct( kt ) ! Dummy routine
	1306	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2848]	1307	WRITE(,) 'dia_dct: You should not have seen this print! error?', kt
	1308	END SUBROUTINE dia_dct
	1309	#endif
	1310
	1311	END MODULE diadct

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

Download in other formats: