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 NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 10761

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

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

Download in other formats: