Context Navigation

diadct.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago
GMED 450 add flush after prints
Property svn:executable set to ``*
File size: 60.1 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 10774

Download in other formats: