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/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 7351

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

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

Download in other formats: