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

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 3281

Last change on this file since 3281 was 3263, checked in by rblod, 12 years ago
Correct cpp name for key_ice_lim see ticket #909
Property svn:executable set to ``*
File size: 53.9 KB

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

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

Download in other formats: