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.

diaharm.F90 in branches/UKMO/r5518_sst_landsea_cpl/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/r5518_sst_landsea_cpl/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 7147

Last change on this file since 7147 was 7147, checked in by jcastill, 8 years ago
Remove svn keywords
File size: 20.0 KB

Rev	Line
[2956]	1	MODULE diaharm
[4292]	2	!!======================================================================
[2956]	3	!! * MODULE diaharm *
	4	!! Harmonic analysis of tidal constituents
[4292]	5	!!======================================================================
	6	!! History : 3.1 ! 2007 (O. Le Galloudec, J. Chanut) Original code
	7	!!----------------------------------------------------------------------
	8	#if defined key_diaharm && defined key_tide
	9	!!----------------------------------------------------------------------
	10	!! 'key_diaharm'
	11	!! 'key_tide'
	12	!!----------------------------------------------------------------------
[2956]	13	USE oce ! ocean dynamics and tracers variables
	14	USE dom_oce ! ocean space and time domain
	15	USE phycst
	16	USE daymod
	17	USE tide_mod
[4683]	18	!
[4292]	19	USE in_out_manager ! I/O units
	20	USE iom ! I/0 library
	21	USE ioipsl ! NetCDF IPSL library
	22	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	23	USE diadimg ! To write dimg
[3294]	24	USE timing ! preformance summary
	25	USE wrk_nemo ! working arrays
[2956]	26
	27	IMPLICIT NONE
	28	PRIVATE
[3294]	29
	30	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .TRUE.
[2956]	31
[3294]	32	INTEGER, PARAMETER :: jpincomax = 2.*jpmax_harmo
	33	INTEGER, PARAMETER :: jpdimsparse = jpincomax30024
[2956]	34
[4683]	35	! !! namelist variables
[4292]	36	INTEGER :: nit000_han ! First time step used for harmonic analysis
	37	INTEGER :: nitend_han ! Last time step used for harmonic analysis
	38	INTEGER :: nstep_han ! Time step frequency for harmonic analysis
[4683]	39	INTEGER :: nb_ana ! Number of harmonics to analyse
[2956]	40
[4292]	41	INTEGER , ALLOCATABLE, DIMENSION(:) :: name
	42	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
	43	REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut , vt , ft
	44	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta , out_u, out_v
[2956]	45
[4292]	46	INTEGER :: ninco, nsparse
	47	INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse
	48	INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1
	49	REAL(wp), DIMENSION(jpdimsparse) :: valuesparse
	50	REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7
	51	REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
	52	REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot
[3294]	53
[4292]	54	CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
[2956]	55
[4292]	56	PUBLIC dia_harm ! routine called by step.F90
[2956]	57
[4292]	58	!!----------------------------------------------------------------------
	59	!! NEMO/OPA 3.5 , NEMO Consortium (2013)
[5215]	60	!! $Id$
[4292]	61	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	62	!!----------------------------------------------------------------------
[2956]	63	CONTAINS
	64
	65	SUBROUTINE dia_harm_init
	66	!!----------------------------------------------------------------------
	67	!! * ROUTINE dia_harm_init *
	68	!!
	69	!! ** Purpose : Initialization of tidal harmonic analysis
	70	!!
[3038]	71	!! ** Method : Initialize frequency array and nodal factor for nit000_han
[2956]	72	!!
	73	!!--------------------------------------------------------------------
	74	INTEGER :: jh, nhan, jk, ji
[4147]	75	INTEGER :: ios ! Local integer output status for namelist read
[2956]	76
[3294]	77	NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
[2956]	78	!!----------------------------------------------------------------------
	79
[3294]	80	IF(lwp) THEN
	81	WRITE(numout,*)
	82	WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
	83	WRITE(numout,*) '~~~~~~~ '
	84	ENDIF
	85	!
	86	CALL tide_init_Wave
	87	!
[4147]	88	REWIND( numnam_ref ) ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
	89	READ ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
	90	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in reference namelist', lwp )
	91
	92	REWIND( numnam_cfg ) ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis
	93	READ ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 )
	94	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist', lwp )
[4624]	95	IF(lwm) WRITE ( numond, nam_diaharm )
[3294]	96	!
	97	IF(lwp) THEN
	98	WRITE(numout,*) 'First time step used for analysis: nit000_han= ', nit000_han
	99	WRITE(numout,*) 'Last time step used for analysis: nitend_han= ', nitend_han
	100	WRITE(numout,*) 'Time step frequency for harmonic analysis: nstep_han= ', nstep_han
[2956]	101	ENDIF
	102
	103	! Basic checks on harmonic analysis time window:
	104	! ----------------------------------------------
[4292]	105	IF( nit000 > nit000_han ) CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000', &
	106	& ' restart capability not implemented' )
	107	IF( nitend < nitend_han ) CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend', &
	108	& 'restart capability not implemented' )
[2956]	109
[4292]	110	IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 ) &
	111	& CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
[2956]	112
[4292]	113	nb_ana = 0
[3294]	114	DO jk=1,jpmax_harmo
	115	DO ji=1,jpmax_harmo
	116	IF(TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN
	117	nb_ana=nb_ana+1
	118	ENDIF
	119	END DO
[4683]	120	END DO
[3294]	121	!
	122	IF(lwp) THEN
	123	WRITE(numout,*) ' Namelist nam_diaharm'
	124	WRITE(numout,*) ' nb_ana = ', nb_ana
	125	CALL flush(numout)
	126	ENDIF
	127	!
	128	IF (nb_ana > jpmax_harmo) THEN
	129	IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nb_ana must be lower than jpmax_harmo, stop'
	130	IF(lwp) WRITE(numout,*) ' jpmax_harmo= ', jpmax_harmo
	131	nstop = nstop + 1
	132	ENDIF
[2956]	133
	134	ALLOCATE(name (nb_ana))
	135	DO jk=1,nb_ana
	136	DO ji=1,jpmax_harmo
	137	IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN
	138	name(jk) = ji
	139	EXIT
	140	END IF
	141	END DO
	142	END DO
	143
	144	! Initialize frequency array:
	145	! ---------------------------
[4292]	146	ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
[2956]	147
[4292]	148	CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
[2956]	149
[2964]	150	IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
[2956]	151
	152	DO jh = 1, nb_ana
	153	IF(lwp) WRITE(numout,*) ' : ',tname(jh),' ',ana_freq(jh)
	154	END DO
	155
	156	! Initialize temporary arrays:
	157	! ----------------------------
[4292]	158	ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
[4683]	159	ana_temp(:,:,:,:) = 0._wp
[2956]	160
	161	END SUBROUTINE dia_harm_init
[4292]	162
	163
[2956]	164	SUBROUTINE dia_harm ( kt )
	165	!!----------------------------------------------------------------------
	166	!! * ROUTINE dia_harm *
	167	!!
	168	!! ** Purpose : Tidal harmonic analysis main routine
	169	!!
[3038]	170	!! ** Action : Sums ssh/u/v over time analysis [nit000_han,nitend_han]
[2956]	171	!!
	172	!!--------------------------------------------------------------------
	173	INTEGER, INTENT( IN ) :: kt
[4292]	174	!
[3294]	175	INTEGER :: ji, jj, jh, jc, nhc
[2956]	176	REAL(wp) :: ztime, ztemp
[3294]	177	!!--------------------------------------------------------------------
	178	IF( nn_timing == 1 ) CALL timing_start('dia_harm')
[2956]	179
[4683]	180	IF( kt == nit000 ) CALL dia_harm_init
[2956]	181
[4683]	182	IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
[2956]	183
[4683]	184	ztime = (kt-nit000+1) * rdt
[2956]	185
[4683]	186	nhc = 0
	187	DO jh = 1, nb_ana
	188	DO jc = 1, 2
	189	nhc = nhc+1
	190	ztemp =( MOD(jc,2) * ft(jh) COS(ana_freq(jh)ztime + vt(jh) + ut(jh)) &
	191	& +(1.-MOD(jc,2))* ft(jh) SIN(ana_freq(jh)ztime + vt(jh) + ut(jh)))
[2956]	192
[4990]	193	DO jj = 1,jpj
	194	DO ji = 1,jpi
[4683]	195	! Elevation
[5586]	196	ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztempsshn(ji,jj)tmask_i(ji,jj)
	197	ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztempun_b(ji,jj)umask_i(ji,jj)
	198	ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztempvn_b(ji,jj)vmask_i(ji,jj)
[4683]	199	END DO
	200	END DO
	201	!
[2956]	202	END DO
[4683]	203	END DO
	204	!
[2956]	205	END IF
	206
[4292]	207	IF ( kt == nitend_han ) CALL dia_harm_end
[2956]	208
[3294]	209	IF( nn_timing == 1 ) CALL timing_stop('dia_harm')
[2956]	210
	211	END SUBROUTINE dia_harm
	212
[4292]	213
[2956]	214	SUBROUTINE dia_harm_end
	215	!!----------------------------------------------------------------------
	216	!! * ROUTINE diaharm_end *
	217	!!
[3038]	218	!! ** Purpose : Compute the Real and Imaginary part of tidal constituents
[2956]	219	!!
[3038]	220	!! ** Action : Decompose the signal on the harmonic constituents
[2956]	221	!!
	222	!!--------------------------------------------------------------------
	223	INTEGER :: ji, jj, jh, jc, jn, nhan, jl
	224	INTEGER :: ksp, kun, keq
	225	REAL(wp) :: ztime, ztime_ini, ztime_end
	226	REAL(wp) :: X1,X2
[3294]	227	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ana_amp
	228	!!--------------------------------------------------------------------
	229	CALL wrk_alloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
[2956]	230
[3294]	231	IF(lwp) WRITE(numout,*)
	232	IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
	233	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
[2956]	234
[3294]	235	ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis
	236	ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis
	237	nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
[2956]	238
[3294]	239	ninco = 2*nb_ana
[2956]	240
[3294]	241	ksp = 0
	242	keq = 0
	243	DO jn = 1, nhan
	244	ztime=( (nhan-jn)ztime_ini + (jn-1)ztime_end )/FLOAT(nhan-1)
	245	keq = keq + 1
	246	kun = 0
[4683]	247	DO jh = 1, nb_ana
	248	DO jc = 1, 2
[3294]	249	kun = kun + 1
	250	ksp = ksp + 1
	251	nisparse(ksp) = keq
	252	njsparse(ksp) = kun
[4292]	253	valuesparse(ksp) = ( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
	254	& + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
[2956]	255	END DO
[3294]	256	END DO
	257	END DO
[2956]	258
[4292]	259	nsparse = ksp
[2956]	260
[3294]	261	! Elevation:
	262	DO jj = 1, jpj
	263	DO ji = 1, jpi
[2956]	264	! Fill input array
[4292]	265	kun = 0
	266	DO jh = 1, nb_ana
	267	DO jc = 1, 2
[3294]	268	kun = kun + 1
	269	ztmp4(kun)=ana_temp(ji,jj,kun,1)
[4292]	270	END DO
	271	END DO
[2956]	272
	273	CALL SUR_DETERMINE(jj)
	274
	275	! Fill output array
	276	DO jh = 1, nb_ana
[3294]	277	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
	278	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
[2956]	279	END DO
[3294]	280	END DO
	281	END DO
[2956]	282
[4292]	283	ALLOCATE( out_eta(jpi,jpj,2*nb_ana), &
	284	& out_u (jpi,jpj,2*nb_ana), &
	285	& out_v (jpi,jpj,2*nb_ana) )
[2956]	286
[3294]	287	DO jj = 1, jpj
	288	DO ji = 1, jpi
[2956]	289	DO jh = 1, nb_ana
[4292]	290	X1 = ana_amp(ji,jj,jh,1)
	291	X2 =-ana_amp(ji,jj,jh,2)
[4990]	292	out_eta(ji,jj,jh ) = X1 * tmask_i(ji,jj)
	293	out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj)
[4683]	294	END DO
	295	END DO
	296	END DO
[2956]	297
[3294]	298	! ubar:
	299	DO jj = 1, jpj
	300	DO ji = 1, jpi
[2956]	301	! Fill input array
	302	kun=0
	303	DO jh = 1,nb_ana
[3294]	304	DO jc = 1,2
	305	kun = kun + 1
	306	ztmp4(kun)=ana_temp(ji,jj,kun,2)
[4683]	307	END DO
	308	END DO
[2956]	309
	310	CALL SUR_DETERMINE(jj+1)
	311
	312	! Fill output array
	313	DO jh = 1, nb_ana
[4683]	314	ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
	315	ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
[2956]	316	END DO
	317
[3294]	318	END DO
	319	END DO
[2956]	320
[3294]	321	DO jj = 1, jpj
	322	DO ji = 1, jpi
[2956]	323	DO jh = 1, nb_ana
[4990]	324	X1= ana_amp(ji,jj,jh,1)
	325	X2=-ana_amp(ji,jj,jh,2)
	326	out_u(ji,jj, jh) = X1 * umask_i(ji,jj)
	327	out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)
	328	ENDDO
	329	ENDDO
	330	ENDDO
[2956]	331
[3294]	332	! vbar:
	333	DO jj = 1, jpj
	334	DO ji = 1, jpi
	335	! Fill input array
	336	kun=0
	337	DO jh = 1,nb_ana
	338	DO jc = 1,2
[2956]	339	kun = kun + 1
[3294]	340	ztmp4(kun)=ana_temp(ji,jj,kun,3)
[4683]	341	END DO
	342	END DO
[2956]	343
[3294]	344	CALL SUR_DETERMINE(jj+1)
[2956]	345
[3294]	346	! Fill output array
	347	DO jh = 1, nb_ana
	348	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
	349	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
	350	END DO
[2956]	351
[3294]	352	END DO
	353	END DO
[2956]	354
[3294]	355	DO jj = 1, jpj
	356	DO ji = 1, jpi
[2956]	357	DO jh = 1, nb_ana
[3294]	358	X1=ana_amp(ji,jj,jh,1)
	359	X2=-ana_amp(ji,jj,jh,2)
[4990]	360	out_v(ji,jj, jh)=X1 * vmask_i(ji,jj)
	361	out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)
[4683]	362	END DO
	363	END DO
	364	END DO
[2956]	365
[3294]	366	CALL dia_wri_harm ! Write results in files
	367	CALL wrk_dealloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
	368	!
	369	END SUBROUTINE dia_harm_end
[2956]	370
[4292]	371
[3294]	372	SUBROUTINE dia_wri_harm
[2956]	373	!!--------------------------------------------------------------------
	374	!! * ROUTINE dia_wri_harm *
	375	!!
	376	!! ** Purpose : Write tidal harmonic analysis results in a netcdf file
	377	!!--------------------------------------------------------------------
	378	CHARACTER(LEN=lc) :: cltext
	379	CHARACTER(LEN=lc) :: &
[3294]	380	cdfile_name_T , & ! name of the file created (T-points)
	381	cdfile_name_U , & ! name of the file created (U-points)
	382	cdfile_name_V ! name of the file created (V-points)
	383	INTEGER :: jh
[2956]	384	!!----------------------------------------------------------------------
	385
	386	#if defined key_dimgout
	387	cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'
	388	cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'
	389	cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'
	390	#endif
	391
	392	IF(lwp) WRITE(numout,*) ' '
[3294]	393	IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
[2956]	394	#if defined key_dimgout
	395	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Output files: ', TRIM(cdfile_name_T)
	396	IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_U)
	397	IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_V)
	398	#endif
	399	IF(lwp) WRITE(numout,*) ' '
	400
	401	! A) Elevation
	402	!/////////////
	403	!
	404	#if defined key_dimgout
	405	cltext='Elevation amplitude and phase'
	406	CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')
	407	#else
	408	DO jh = 1, nb_ana
	409	CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
	410	CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
	411	END DO
	412	#endif
	413
	414	! B) ubar
	415	!/////////
	416	!
	417	#if defined key_dimgout
	418	cltext='ubar amplitude and phase'
	419	CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')
	420	#else
	421	DO jh = 1, nb_ana
	422	CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
	423	CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) )
	424	END DO
	425	#endif
	426
	427	! C) vbar
	428	!/////////
	429	!
	430	#if defined key_dimgout
	431	cltext='vbar amplitude and phase'
	432	CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')
	433	#else
	434	DO jh = 1, nb_ana
[4683]	435	CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh ) )
	436	CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
[2956]	437	END DO
	438	#endif
[4683]	439	!
[3294]	440	END SUBROUTINE dia_wri_harm
[2956]	441
[4292]	442
[3294]	443	SUBROUTINE SUR_DETERMINE(init)
[4683]	444	!!---------------------------------------------------------------------------------
	445	!! * ROUTINE SUR_DETERMINE *
	446	!!
	447	!!
	448	!!
	449	!!---------------------------------------------------------------------------------
	450	INTEGER, INTENT(in) :: init
	451	!
	452	INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
	453	REAL(wp) :: zval1, zval2, zx1
	454	REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
	455	INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
	456	!---------------------------------------------------------------------------------
	457	CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
	458	CALL wrk_alloc( jpincomax , ipos2 , ipivot )
[3294]	459
[4683]	460	IF( init == 1 ) THEN
	461	IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
	462	IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
	463	!
	464	ztmp3(:,:) = 0._wp
	465	!
	466	DO jk1_sd = 1, nsparse
	467	DO jk2_sd = 1, nsparse
	468	nisparse(jk2_sd) = nisparse(jk2_sd)
	469	njsparse(jk2_sd) = njsparse(jk2_sd)
	470	IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
	471	ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) &
	472	& + valuesparse(jk1_sd)*valuesparse(jk2_sd)
	473	ENDIF
	474	END DO
[4292]	475	END DO
[4683]	476	!
	477	DO jj_sd = 1 ,ninco
	478	ipos1(jj_sd) = jj_sd
	479	ipos2(jj_sd) = jj_sd
	480	END DO
	481	!
	482	DO ji_sd = 1 , ninco
	483	!
	484	!find greatest non-zero pivot:
	485	zval1 = ABS(ztmp3(ji_sd,ji_sd))
	486	!
	487	ipivot(ji_sd) = ji_sd
	488	DO jj_sd = ji_sd, ninco
	489	zval2 = ABS(ztmp3(ji_sd,jj_sd))
	490	IF( zval2.GE.zval1 )THEN
	491	ipivot(ji_sd) = jj_sd
	492	zval1 = zval2
	493	ENDIF
	494	END DO
	495	!
	496	DO ji1_sd = 1, ninco
	497	zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd)
	498	zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd))
	499	ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd)
	500	ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
	501	END DO
	502	!
	503	ipos2(ji_sd) = ipos1(ipivot(ji_sd))
	504	ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
	505	ipos1(ji_sd) = ipos2(ji_sd)
	506	ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
	507	zpivot(ji_sd) = ztmp3(ji_sd,ji_sd)
	508	DO jj_sd = 1, ninco
	509	ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
	510	END DO
	511	!
	512	DO ji2_sd = ji_sd+1, ninco
	513	zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
	514	DO jj_sd=1,ninco
	515	ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
	516	END DO
	517	END DO
	518	!
	519	END DO
	520	!
	521	ENDIF ! End init==1
[3294]	522
[4683]	523	DO ji_sd = 1, ninco
	524	ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
[3294]	525	DO ji2_sd = ji_sd+1, ninco
[4683]	526	ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
	527	END DO
	528	END DO
[3294]	529
[4683]	530	!system solving:
	531	ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
	532	ji_sd = ninco
	533	DO ji_sd = ninco-1, 1, -1
	534	zx1 = 0._wp
	535	DO jj_sd = ji_sd+1, ninco
	536	zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
	537	END DO
	538	ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
	539	END DO
[3294]	540
[4683]	541	DO jj_sd =1, ninco
	542	ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
	543	END DO
[3294]	544
[4683]	545	CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
	546	CALL wrk_dealloc( jpincomax , ipos2 , ipivot )
	547	!
	548	END SUBROUTINE SUR_DETERMINE
[3294]	549
[2956]	550	#else
	551	!!----------------------------------------------------------------------
	552	!! Default case : Empty module
	553	!!----------------------------------------------------------------------
	554	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .FALSE.
	555	CONTAINS
	556	SUBROUTINE dia_harm ( kt ) ! Empty routine
	557	INTEGER, INTENT( IN ) :: kt
	558	WRITE(,) 'dia_harm: you should not have seen this print'
	559	END SUBROUTINE dia_harm
[4292]	560	#endif
[2956]	561
	562	!!======================================================================
	563	END MODULE diaharm

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

Download in other formats: