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/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DIA/diaharm.F90 @ 10826

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

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

Download in other formats: