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

Context Navigation

source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 11132

Last change on this file since 11132 was 11132, checked in by jcastill, 5 years ago
Remove svn keywords
Property svn:executable set to ``*
File size: 59.9 KB

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

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

Download in other formats: