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

Last change on this file since 2864 was 2864, checked in by cbricaud, 13 years ago

add dia_dct_init Dummy routine

Property svn:executable set to *

File size: 47.0 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
[2864]	337	!WRITE(narea+200,*)' # I J : ',iiloc,ijloc
[2848]	338	ENDIF
	339
	340	ENDDO
	341
	342	secs(jsec)%nb_point=iptloc !store number of section's points
	343
	344	!debug
	345	!-----
[2864]	346	!IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
	347	IF( ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
[2848]	348	WRITE(numout,*)" List of points selected by the proc:"
	349	DO jpt = 1,iptloc
	350	iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
	351	ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
	352	WRITE(numout,*)' # I J : ',iiglo,ijglo
[2864]	353	!WRITE(narea+200,*)' # I J : ',iiglo,ijglo
[2848]	354	ENDDO
	355	ENDIF
	356
	357	!remove redundant points between processors
	358	!look NEMO documentation for more explanations
	359	! lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
	360	lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) ) lldebug = .TRUE.
	361	IF( iptloc .NE. 0 )THEN
	362	CALL removepoints(secs(jsec),'I','top_list',lldebug)
	363	CALL removepoints(secs(jsec),'I','bot_list',lldebug)
	364	CALL removepoints(secs(jsec),'J','top_list',lldebug)
	365	CALL removepoints(secs(jsec),'J','bot_list',lldebug)
	366	ENDIF
	367
	368	!debug
	369	!-----
	370	IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
	371	WRITE(numout,*)" List of points after removepoints:"
	372	DO jpt = 1,iptloc
	373	iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
	374	ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
	375	WRITE(numout,*)' # I J : ',iiglo,ijglo
	376	ENDDO
	377	ENDIF
	378
	379	ELSE ! iptglo = 0
	380	IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )&
	381	WRITE(numout,*)' Not points for this section.'
	382	ENDIF
	383
	384	ENDDO !end of the loop on jsec
	385
	386	nb_sec = jsec-1 !number of section read in the file
	387
	388	END SUBROUTINE readsec
	389
	390	SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
	391	!!---------------------------------------------------------------------------
	392	!! *** function removepoints
	393	!!
[2854]	394	!! ** Purpose ::
	395	!! remove points which are common to 2 procs
	396	!!
	397	!!
[2848]	398	!----------------------------------------------------------------------------
	399	!! * arguments
	400	TYPE(SECTION),INTENT(INOUT) :: sec
	401	CHARACTER(len=1),INTENT(IN) :: cdind ! = 'I'/'J'
	402	CHARACTER(len=8),INTENT(IN) :: cdextr ! = 'top_list'/'bot_list'
	403	LOGICAL,INTENT(IN) :: ld_debug
	404
	405	!! * Local variables
	406	INTEGER :: iextr ,& !extremity of listpoint that we verify
	407	iind ,& !coord of listpoint that we verify
	408	itest ,& !indice value of the side of the domain
	409	!where points could be redundant
	410	isgn ,& !way of course in listpoint
	411	istart,iend !first and last points selected in listpoint
	412	INTEGER :: jpoint =0 !loop on list points
	413	INTEGER,DIMENSION(nb_point_max) :: idirec !contains temporare sec%direction
	414	INTEGER,DIMENSION(2,nb_point_max) :: icoord !contains temporare sec%listpoint
	415	!----------------------------------------------------------------------------
	416	IF( ld_debug )WRITE(numout,*)' -------------------------'
	417	IF( ld_debug )WRITE(numout,*)' removepoints in listpoint'
	418
	419	!iextr=extremity of list_point that we verify
	420	IF ( cdextr=='bot_list' )THEN ; iextr=1 ; isgn=1
	421	ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
	422	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdextr")
	423	ENDIF
	424
	425	!which coordinate shall we verify ?
	426	IF ( cdind=='I' )THEN ; itest=nlei ; iind=1
	427	ELSE IF ( cdind=='J' )THEN ; itest=nlej ; iind=2
	428	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdind")
	429	ENDIF
	430
	431	IF( ld_debug )THEN
	432	WRITE(numout,*)' case: coord/list extr/domain side'
	433	WRITE(numout,*)' ', cdind,' ',cdextr,' ',itest
	434	WRITE(numout,*)' Actual number of points: ',sec%nb_point
	435	ENDIF
	436
	437	icoord(1,1:nb_point_max) = sec%listPoint%I
	438	icoord(2,1:nb_point_max) = sec%listPoint%J
	439	idirec = sec%direction
	440	sec%listPoint = POINT_SECTION(0,0)
	441	sec%direction = 0
	442
	443
	444	jpoint=iextr+isgn
	445	DO WHILE ( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point .AND. &
	446	icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest)
	447	jpoint=jpoint+isgn
	448	ENDDO
	449	IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
	450	ELSE ; istart=1 ; iend=jpoint+1
	451	ENDIF
	452	sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
	453	sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
	454	sec%direction(1:1+iend-istart) = idirec(istart:iend)
	455	sec%nb_point = iend-istart+1
	456
	457	IF( ld_debug )THEN
	458	WRITE(numout,*)' Number of points after removepoints :',sec%nb_point
	459	WRITE(numout,*)' sec%direction after removepoints :',sec%direction(1:sec%nb_point)
	460	ENDIF
	461
	462	END SUBROUTINE removepoints
	463
	464	SUBROUTINE transport(sec,ld_debug)
	465	!!---------------------------------------------------------------------
	466	!! * ROUTINE transport *
	467	!!
	468	!! ** Purpose : Compute the transport trough a sec
	469	!!
	470	!! ** Method :Transport through a given section is equal to the sum of transports
	471	!! computed on each proc.
	472	!! On each proc,transport is equal to the sum of transport computed through
	473	!! segments linking each points of sec%listPoint with the next one.
	474	!! There are several loops:
	475	!! loop on the density/temperature/salinity/level classes
	476	!! loop on the segment between 2 nodes
	477	!! loop on the level jk
	478	!! test on the density/temperature/salinity/level
	479	!!
	480	!! ** Output: sec%transport: transport in the 2 direction and temperature,
	481	!! salinity, density meaned by the transport
[2854]	482	!!
	483	!!
[2848]	484	!!-------------------------------------------------------------------
	485	!! * Arguments
	486	TYPE(SECTION),INTENT(INOUT) :: sec
	487	LOGICAL ,INTENT(IN) :: ld_debug
	488
	489	!! * Local variables
	490	INTEGER :: jk,jseg,jclass, &!loop on level/segment/classes
	491	iptegrid, &!=0 if section is constant along i/j
	492	isgnu , isgnv !
	493	INTEGER :: ii, ij ! local integer
	494	REAL(wp):: zumid , zvmid ,&!U/V velocity on a cell segment
	495	zumid_ice , zvmid_ice ,&!U/V ice velocity
	496	zTnorm ,&!transport of velocity through one cell's sides
	497	ztransp1 , ztransp2 ,&!total transport in sens 1 and 2
	498	ztemp1 , ztemp2 ,&!temperature transport "
	499	zrhoi1 , zrhoi2 ,&!mass transport "
	500	zrhop1 , zrhop2 ,&!mass transport "
	501	zsal1 , zsal2 ,&!salinity transport "
	502	zice_vol_pos , zice_vol_neg ,&!volumic ice transport "
	503	zice_surf_pos, zice_surf_neg !surfacic ice transport "
	504	REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
	505	REAL(wp) :: aj,i0,j0,i1,j1,i,j
	506
	507	TYPE(POINT_SECTION) :: k
	508	INTEGER,DIMENSION(1):: ish
	509	INTEGER,DIMENSION(2):: ish2
	510	REAL(wp),DIMENSION(nb_type_class,nb_class_max)::zsum
	511	REAL(wp),DIMENSION(nb_type_class*nb_class_max)::zwork ! temporary vector for mpp_sum
	512	!!--------------------------------------------------------
	513
	514	IF( ld_debug )WRITE(numout,*)' Compute transport'
	515
	516	!----------------!
	517	! INITIALIZATION !
	518	!----------------!
	519	zsum = 0._wp
	520	zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp
	521	zice_vol_pos = 0._wp ; zice_vol_neg = 0._wp
	522
	523	!compute iptegrid: iptegrid=0 if section is constant along i or j
	524	iptegrid = 0
	525	jseg = 1
	526	DO WHILE( (jseg .LT. MAX(sec%nb_point-1,0)) .AND. iptegrid .NE. 1 )
	527	IF(sec%direction(jseg) .NE. sec%direction(1)) iptegrid = 1
	528	jseg = jseg+1
	529	ENDDO
	530	IF(lk_mpp) CALL mpp_sum(iptegrid)
	531
	532	!---------------------------!
	533	! COMPUTE TRANSPORT !
	534	!---------------------------!
	535	IF(sec%nb_point .NE. 0)THEN
	536
	537	!----------------------------------------------------------------------------------------------------
	538	!Compute sign for velocities:
	539	!
	540	!convention:
	541	! non horizontal section: sens + is toward left hand of section
	542	! horizontal section: sens + is toward north of section
	543	!
	544	!
	545	! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0
	546	! ---------------- ----------------- --------------- --------------
	547	!
	548	! isgnv=1 sens +
	549	! ______ _____ ______(i1,j1) (i1,j1)
	550	!(i0,j0) \| //\| \| \| sens +
	551	! \| isgnu=1 // \| \|isgnu=1 \|isgnu=1 /\|\
	552	! \|_______ // ______\| \\ \| ---\ \|
	553	! \| \| isgnv=-1 \\ \| \| ---/ sens + ____________(i1,j1)
	554	! \| \| __\\\| \| (i0,j0)
	555	! \|(i1,j1) \|(i0,j0) sens + \| isgnv=1
	556	! (i0,j0)
	557	!----------------------------------------------------------------------------------------------------
	558	isgnu = 1
	559	IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1
	560	ELSE ; isgnv = 1
	561	ENDIF
	562
	563	IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
	564
	565	!--------------------------------------!
	566	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
	567	!--------------------------------------!
	568	DO jseg=1,MAX(sec%nb_point-1,0)
	569
	570	!------------------------------------------
	571	! Select good coordinate to have velocity of the segment
	572	!
	573	! \|U CASE(3)
	574	! CASE(0) \| j+1
	575	! ____V_____\|F____V_CASE(1)__
	576	! \|
	577	! \| j
	578	! T \|U CASE(2)
	579	! i \| i+1
	580	!-------------------------------------------
	581	SELECT CASE( sec%direction(jseg) )
	582	CASE(0) ; k = sec%listPoint(jseg)
	583	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
	584	CASE(2) ; k = sec%listPoint(jseg)
	585	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
	586	END SELECT
[2864]	587
[2848]	588	!-------------------------------
	589	! LOOP ON THE DENSITY CLASSES \|
	590	!-------------------------------
	591	!The computation is made for each density class
	592	DO jclass=1,MAX(1,sec%nb_class-1)
	593
	594	ztransp1=0._wp ; zrhoi1=0._wp ; zrhop1=0._wp ; ztemp1=0._wp ;zsal1=0._wp
	595	ztransp2=0._wp ; zrhoi2=0._wp ; zrhop2=0._wp ; ztemp2=0._wp ;zsal2=0._wp
	596
	597	!---------------------------\|
	598	! LOOP ON THE LEVEL \|
	599	!---------------------------\|
	600	!Sum of the transport on the vertical
	601	DO jk=1,jpk
	602
	603
	604	! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point
	605	SELECT CASE( sec%direction(jseg) )
	606	CASE(0,1)
	607	ztn = interp(k%I,k%J,jk,'V',tn)
	608	zsn = interp(k%I,k%J,jk,'V',sn)
	609	zrhop = interp(k%I,k%J,jk,'V',rhop)
	610	zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
	611	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) vmask(k%I,k%J,1)
	612	CASE(2,3)
	613	ztn = interp(k%I,k%J,jk,'U',tn)
	614	zsn = interp(k%I,k%J,jk,'U',sn)
	615	zrhop = interp(k%I,k%J,jk,'U',rhop)
	616	zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
	617	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) umask(k%I,k%J,1)
	618	END SELECT
	619
	620	zfsdep= gdept(k%I,k%J,jk)
	621
	622	!----------------------------------------------!
	623	!TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
	624	!----------------------------------------------!
[2864]	625	IF(ld_debug)write(narea+200,*)k ,jk
[2848]	626
	627	IF ( ( ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. &
	628	( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. &
	629	( sec%zsigp(jclass) .EQ. 99.)) .AND. &
[2864]	630	((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. &
	631	( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. &
[2848]	632	( sec%zsigi(jclass) .EQ. 99.)) .AND. &
	633	((( zsn .GT. sec%zsal(jclass)) .AND. &
	634	( zsn .LE. sec%zsal(jclass+1))) .OR. &
	635	( sec%zsal(jclass) .EQ. 99.)) .AND. &
	636	((( ztn .GE. sec%ztem(jclass)) .AND. &
	637	( ztn .LE. sec%ztem(jclass+1))) .OR. &
	638	( sec%ztem(jclass) .EQ.99.)) .AND. &
	639	((( zfsdep .GE. sec%zlay(jclass)) .AND. &
	640	( zfsdep .LE. sec%zlay(jclass+1))) .OR. &
	641	( sec%zlay(jclass) .EQ. 99. )))) THEN
	642
	643
	644	!compute velocity with the good sens
	645	SELECT CASE( sec%direction(jseg) )
	646	CASE(0,1)
	647	zumid=0.
	648	zvmid=isgnvvn(k%I,k%J,jk)vmask(k%I,k%J,jk)
	649	CASE(2,3)
	650	zumid=isgnuun(k%I,k%J,jk)umask(k%I,k%J,jk)
	651	zvmid=0.
	652	END SELECT
	653
	654	!velocity* cell's lenght * cell's thickness
	655	zTnorm=zumide2u(k%I,k%J) fse3u(k%I,k%J,jk)+ &
	656	zvmide1v(k%I,k%J) fse3v(k%I,k%J,jk)
	657
	658	#if ! defined key_vvl
	659	!add transport due to free surface
	660	IF( jk==1 )THEN
	661	zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
	662	zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
	663	ENDIF
	664	#endif
[2864]	665	!IF(ld_debug)write(narea+200,*)k ,jk,zTnorm
[2848]	666	!COMPUTE TRANSPORT
	667	!zTnorm=transport through one cell for one class
	668	!ztransp1 or ztransp2=transport through one cell i
	669	! for one class for one direction
	670	IF( zTnorm .GE. 0 )THEN
	671
	672	ztransp1=zTnorm+ztransp1
	673
	674	IF ( sec%llstrpond ) THEN
	675	ztemp1 = ztemp1 + zTnorm * ztn
	676	zsal1 = zsal1 + zTnorm * zsn
	677	zrhoi1 = zrhoi1 + zTnorm * zrhoi
	678	zrhop1 = zrhop1 + zTnorm * zrhop
	679	ENDIF
	680
	681	ELSE
	682
	683	ztransp2=(zTnorm)+ztransp2
	684
	685	IF ( sec%llstrpond ) THEN
	686	ztemp2 = ztemp2 + zTnorm * ztn
	687	zsal2 = zsal2 + zTnorm * zsn
	688	zrhoi2 = zrhoi2 + zTnorm * zrhoi
	689	zrhop2 = zrhop2 + zTnorm * zrhop
	690	ENDIF
	691	ENDIF
	692
	693
	694	ENDIF ! end of density test
	695	ENDDO!end of loop on the level
	696
	697	!ZSUM=TRANSPORT FOR EACH CLASSES FOR THE DIRECTIONS
	698	!---------------------------------------------------
	699	zsum(1,jclass) = zsum(1,jclass)+ztransp1
	700	zsum(2,jclass) = zsum(2,jclass)+ztransp2
	701	IF( sec%llstrpond )THEN
	702	zsum(3 ,jclass) = zsum( 3,jclass)+zrhoi1
	703	zsum(4 ,jclass) = zsum( 4,jclass)+zrhoi2
	704	zsum(5 ,jclass) = zsum( 5,jclass)+zrhop1
	705	zsum(6 ,jclass) = zsum( 6,jclass)+zrhop2
	706	zsum(7 ,jclass) = zsum( 7,jclass)+ztemp1
	707	zsum(8 ,jclass) = zsum( 8,jclass)+ztemp2
	708	zsum(9 ,jclass) = zsum( 9,jclass)+zsal1
	709	zsum(10,jclass) = zsum(10,jclass)+zsal2
	710	ENDIF
	711
	712	ENDDO !end of loop on the density classes
	713
	714	#if defined key_ice_lim
	715
	716	!ICE CASE
	717	!------------
	718	IF( sec%ll_ice_section )THEN
	719	SELECT CASE (sec%direction(jseg))
	720	CASE(0)
	721	zumid_ice = 0
	722	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
	723	CASE(1)
	724	zumid_ice = 0
	725	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
	726	CASE(2)
	727	zvmid_ice = 0
	728	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
	729	CASE(3)
	730	zvmid_ice = 0
	731	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
	732	END SELECT
	733
	734	zTnorm=zumid_icee2u(k%I,k%J)+zvmid_icee1v(k%I,k%J)
	735
	736	IF( zTnorm .GE. 0)THEN
	737	zice_vol_pos = (zTnorm)* &
	738	(1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
	739	*(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + &
	740	hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
	741	+zice_vol_pos
	742	zice_surf_pos = (zTnorm)* &
	743	(1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
	744	+zice_surf_pos
	745	ELSE
	746	zice_vol_neg=(zTnorm)* &
	747	(1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
	748	*(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + &
	749	hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
	750	+zice_vol_neg
	751	zice_surf_neg=(zTnorm)* &
	752	(1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
	753	+zice_surf_neg
	754	ENDIF
	755
	756	zsum(11,1) = zsum(11,1)+zice_vol_pos
	757	zsum(12,1) = zsum(12,1)+zice_vol_neg
	758	zsum(13,1) = zsum(13,1)+zice_surf_pos
	759	zsum(14,1) = zsum(14,1)+zice_surf_neg
	760
	761	ENDIF !end of ice case
	762	#endif
	763
	764	ENDDO !end of loop on the segment
	765
	766
	767	ELSE !if nb_inmesh=0
	768	zsum(1:2,:)=0.
	769	IF (sec%llstrpond) zsum(3:10,:)=0.
	770	zsum( 11:14,:)=0.
	771	ENDIF !end of nb_inmesh=0 case
	772
	773	!------------------!
	774	!SUM ON ALL PROCS !
	775	!------------------!
	776	IF( lk_mpp )THEN
	777	ish(1) = nb_type_class*nb_class_max ; ish2(1)=nb_type_class ; ish2(2)=nb_class_max
	778	zwork(:)= RESHAPE(zsum(:,:), ish )
	779	CALL mpp_sum(zwork, ish(1))
	780	zsum(:,:)= RESHAPE(zwork,ish2)
	781	ENDIF
	782
	783	!-------------------------------\|
	784	!FINISH COMPUTING TRANSPORTS \|
	785	!-------------------------------\|
	786	DO jclass=1,MAX(1,sec%nb_class-1)
	787	sec%transport(1,jclass)=sec%transport(1,jclass)+zsum(1,jclass)*1.E-6
	788	sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6
	789	IF( sec%llstrpond ) THEN
	790	IF( zsum(1,jclass) .NE. 0 ) THEN
	791	sec%transport(3,jclass)=sec%transport(3,jclass)+zsum(3,jclass)/zsum(1,jclass)
	792	sec%transport(5,jclass)=sec%transport(5,jclass)+zsum(5,jclass)/zsum(1,jclass)
[2854]	793	sec%transport(7,jclass)=sec%transport(7,jclass)+zsum(7,jclass)
	794	sec%transport(9,jclass)=sec%transport(9,jclass)+zsum(9,jclass)
[2848]	795	ELSE
	796	sec%transport(3,jclass)=0.
	797	sec%transport(5,jclass)=0.
	798	sec%transport(7,jclass)=0.
	799	sec%transport(9,jclass)=0.
	800	ENDIF
	801	IF( zsum(2,jclass) .NE. 0 )THEN
	802	sec%transport( 4,jclass)=sec%transport( 4,jclass)+zsum( 4,jclass)/zsum(2,jclass)
	803	sec%transport( 6,jclass)=sec%transport( 6,jclass)+zsum( 6,jclass)/zsum(2,jclass)
[2854]	804	sec%transport( 8,jclass)=sec%transport( 8,jclass)+zsum( 8,jclass)
	805	sec%transport(10,jclass)=sec%transport(10,jclass)+zsum(10,jclass)
[2848]	806	ELSE
	807	sec%transport( 4,jclass)=0.
	808	sec%transport( 6,jclass)=0.
	809	sec%transport( 8,jclass)=0.
	810	sec%transport(10,jclass)=0.
	811	ENDIF
	812	ELSE
	813	sec%transport( 3,jclass)=0.
	814	sec%transport( 4,jclass)=0.
	815	sec%transport( 5,jclass)=0.
	816	sec%transport( 6,jclass)=0.
	817	sec%transport( 7,jclass)=0.
	818	sec%transport( 8,jclass)=0.
	819	sec%transport(10,jclass)=0.
	820	ENDIF
	821	ENDDO
	822
	823	IF( sec%ll_ice_section ) THEN
	824	sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6
	825	sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6
	826	sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6
	827	sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6
	828	ENDIF
	829
	830	END SUBROUTINE transport
	831
	832	SUBROUTINE dia_dct_wri(kt,jsec,sec)
	833	!!-------------------------------------------------------------
	834	!! Write transport output in numdct
	835	!!
[2854]	836	!! Purpose: Write transports in ascii files
	837	!!
	838	!! Method:
	839	!! 1. Write volume transports in "volume_transport"
	840	!! Unit: Sv : Surface * Velocity / 1.e6
	841	!!
	842	!! 2. Write heat transports in "heat_transport"
	843	!! Unit: Peta W : Surface * Velocity * T * rhau * Cp / 1.e15
	844	!!
	845	!! 3. Write salt transports in "salt_transport"
	846	!! Unit: Mega m^3 / s : Surface * Velocity * S / 1.e6
	847	!!
[2848]	848	!!-------------------------------------------------------------
	849	!!arguments
	850	INTEGER, INTENT(IN) :: kt
	851	TYPE(SECTION), INTENT(INOUT) :: sec
	852	INTEGER ,INTENT(IN) :: jsec
	853
	854	!!local declarations
[2854]	855	REAL(wp) ,DIMENSION(nb_type_class):: zsumclass
[2848]	856	INTEGER :: jcl,ji !Dummy loop
	857	CHARACTER(len=2) :: classe
	858	REAL(wp) :: zbnd1,zbnd2
	859	!!-------------------------------------------------------------
	860
	861	zsumclass(:)=0._wp
	862
	863	DO jcl=1,MAX(1,sec%nb_class-1)
	864
	865	! Mean computation
	866	sec%transport(:,jcl)=sec%transport(:,jcl)/(nn_dctwri/nn_dct)
	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)1000._wprcp/1.e15,sec%transport(8,jcl)1000._wprcp/1.e15, &
	921	( sec%transport(7,jcl)+sec%transport(8,jcl) )1000._wprcp/1.e15
[2848]	922	!write salt transport per class
	923	WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
	924	jcl,classe,zbnd1,zbnd2,&
[2857]	925	sec%transport(9,jcl)1000._wp/1.e9,sec%transport(10,jcl)1000._wp/1.e9,&
	926	(sec%transport(9,jcl)+sec%transport(10,jcl))*1000._wp/1.e9
[2848]	927	ENDIF
	928
	929	ENDDO
	930
	931	zbnd1 = 0._wp
	932	zbnd2 = 0._wp
	933	jcl=0
	934
	935	!write total volume transport
	936	WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection, &
	937	jcl,"total",zbnd1,zbnd2,&
	938	zsumclass(1),zsumclass(2),zsumclass(1)+zsumclass(2)
	939
	940	IF( sec%llstrpond )THEN
	941
	942	!write total heat transport
	943	WRITE(numdct_heat,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
	944	jcl,"total",zbnd1,zbnd2,&
[2854]	945	zsumclass(7)* 1000._wprcp/1.e15,zsumclass(8) 1000._wprcp/1.e15,(zsumclass(7)+zsumclass(8) ) 1000._wp*rcp/1.e15
[2848]	946	!write total salt transport
	947	WRITE(numdct_salt,119) ndastp,kt,jsec,sec%name,sec%slopeSection, &
	948	jcl,"total",zbnd1,zbnd2,&
[2857]	949	zsumclass(9)1000._wp/1.e9,zsumclass(10)1000._wp/1.e9,(zsumclass(9)+zsumclass(10))*1000._wp/1.e9
[2848]	950	ENDIF
	951
	952
	953	IF ( sec%ll_ice_section) THEN
	954	!write total ice volume transport
	955	WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,&
	956	jcl,"ice_vol",zbnd1,zbnd2,&
	957	sec%transport(9,1),sec%transport(10,1),&
	958	sec%transport(9,1)+sec%transport(10,1)
	959	!write total ice surface transport
	960	WRITE(numdct_vol,118) ndastp,kt,jsec,sec%name,sec%slopeSection,&
	961	jcl,"ice_surf",zbnd1,zbnd2,&
	962	sec%transport(11,1),sec%transport(12,1), &
	963	sec%transport(11,1)+sec%transport(12,1)
	964	ENDIF
	965
	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	END SUBROUTINE dia_dct_wri
	970
	971	FUNCTION interp(ki, kj, kk, cd_point, ptab)
	972	!!----------------------------------------------------------------------
	973	!!
[2854]	974	!! Purpose: compute Temperature/Salinity/density on U-point or V-point
[2848]	975	!!
[2854]	976	!! \| I \| I+1 \|
	977	!! \| \| \|
	978	!! -------------------------------- 1.Compute Tbis:
	979	!! \| \| \| interpolate T(I,J,K) et T(I,J,K+1)
	980	!! \| \| \|
	981	!! \| \| \| 2. 1.Compute Temp/Salt/Rho à U point:
	982	!! K \| T \| \| interpolate Tbis et T(I+U,J,K+1)
	983	!! \| \| \|
	984	!! \| \| \|
	985	!! \| \| \|
	986	!! --------------------------------
	987	!! \| \| \|
	988	!! \| \| \|
	989	!! \| \| \|
	990	!!K+1 \| Tbis U T \|
	991	!! \| \| \|
	992	!! \| T \| \|
	993	!! \| \|------------\|
	994	!! \| \| partials \|
	995	!! \| \| steps \|
	996	!! --------------------------------
[2848]	997	!!
	998	!!----------------------------------------------------------------------
	999	!*arguments
	1000	INTEGER, INTENT(IN) :: ki, kj, kk
	1001	CHARACTER(len=1), INTENT(IN) :: cd_point
	1002	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab
	1003	REAL(wp) :: interp
	1004
	1005	!*local declations
	1006	INTEGER :: ii1, ij1, ii2, ij2
	1007	REAL(wp):: ze3w, zfse3, zmax1, zmax2, zbis
	1008	!!----------------------------------------------------------------------
	1009
	1010	IF( cd_point=='U' )THEN
	1011	ii1 = ki ; ij1 = kj
	1012	ii2 = ki+1 ; ij2 = kj
	1013	ELSE ! cd_point=='V'
	1014	ii1 = ki ; ij1 = kj
	1015	ii2 = ki ; ij2 = kj+1
	1016	ENDIF
	1017
	1018	#if defined key_vvl
	1019
	1020	ze3w = fsve3w(ii2,ij2,kk) - fsve3w(ii1,ij1,kk)
	1021
	1022	zfse3 = fsve3w(ii1,ij1,kk) * ( 1 + sshn(ii2,ij2) * mut(ii2,ij2,kk) )
	1023	zmax1 = ( fse3w(ii2,ij2,kk) - zfse3 )/ fse3w(ii2,ij2,kk)
	1024
	1025	zfse3 = fsve3w(ii2,ij2,kk) * ( 1 + sshn(ii1,ij1) * mut(ii1,ij1,kk) )
	1026	zmax2 = ( fse3w(ii1,ij1,kk) - zfse3 )/ fse3w(ii1,ij1,kk)
	1027
	1028	#else
	1029
	1030	ze3w = fse3w(ii2,ij2,kk)-fse3w(ii1,ij1,kk)
	1031	zmax1 = ze3w / fse3w(ii2,ij2,kk)
	1032	zmax2 = -ze3w / fse3w(ii1,ij1,kk)
	1033
	1034	#endif
	1035
	1036	IF(kk .NE. 1)THEN
	1037
	1038	IF( ze3w >= 0. )THEN
	1039	!zbis
	1040	zbis = ptab(ii2,ij2,kk) + zmax1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )
	1041	! result
	1042	interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + zbis )
	1043	ELSE
	1044	!zbis
	1045	zbis = ptab(ii1,ij1,kk) + zmax2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
	1046	! result
	1047	interp = umask(ii1,ij1,kk) * 0.5*( zbis + ptab(ii2,ij2,kk) )
	1048	ENDIF
	1049
	1050	ELSE
	1051	interp = umask(ii1,ij1,kk) * 0.5*( ptab(ii1,ij1,kk) + ptab(ii2,ij2,kk) )
	1052	ENDIF
	1053
	1054	END FUNCTION interp
	1055
	1056	#else
	1057	!!----------------------------------------------------------------------
	1058	!! Default option : Dummy module
	1059	!!----------------------------------------------------------------------
	1060	LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .FALSE. !: diamht flag
[2864]	1061	PUBLIC
[2848]	1062	CONTAINS
[2864]	1063
	1064	SUBROUTINE dia_dct_init ! Dummy routine
	1065	WRITE(,) 'dia_dct_init: You should not have seen this print! error?', kt
	1066	END SUBROUTINE dia_dct_init
	1067
[2848]	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: