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

Context Navigation

source: branches/2011/dev_r2802_MERCATOR10_diadct/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 2854

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

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

Download in other formats: