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 @ 2848

Last change on this file since 2848 was 2848, checked in by cbricaud, 13 years ago
add routines for online transports computing
Property svn:executable set to ``*
File size: 45.5 KB

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

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

Download in other formats: