Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

diaharm.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 18.7 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
[2956]	24
	25	IMPLICIT NONE
	26	PRIVATE
[3294]	27
	28	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .TRUE.
[2956]	29
[3294]	30	INTEGER, PARAMETER :: jpincomax = 2.*jpmax_harmo
	31	INTEGER, PARAMETER :: jpdimsparse = jpincomax30024
[2956]	32
[4683]	33	! !! namelist variables
[4292]	34	INTEGER :: nit000_han ! First time step used for harmonic analysis
	35	INTEGER :: nitend_han ! Last time step used for harmonic analysis
	36	INTEGER :: nstep_han ! Time step frequency for harmonic analysis
[4683]	37	INTEGER :: nb_ana ! Number of harmonics to analyse
[2956]	38
[4292]	39	INTEGER , ALLOCATABLE, DIMENSION(:) :: name
	40	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
	41	REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut , vt , ft
	42	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta , out_u, out_v
[2956]	43
[4292]	44	INTEGER :: ninco, nsparse
	45	INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse
	46	INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1
	47	REAL(wp), DIMENSION(jpdimsparse) :: valuesparse
	48	REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7
	49	REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
	50	REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot
[3294]	51
[4292]	52	CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
[2956]	53
[4292]	54	PUBLIC dia_harm ! routine called by step.F90
[2956]	55
[4292]	56	!!----------------------------------------------------------------------
	57	!! NEMO/OPA 3.5 , NEMO Consortium (2013)
[5215]	58	!! $Id$
[4292]	59	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	60	!!----------------------------------------------------------------------
[2956]	61	CONTAINS
	62
	63	SUBROUTINE dia_harm_init
	64	!!----------------------------------------------------------------------
	65	!! * ROUTINE dia_harm_init *
	66	!!
	67	!! ** Purpose : Initialization of tidal harmonic analysis
	68	!!
[3038]	69	!! ** Method : Initialize frequency array and nodal factor for nit000_han
[2956]	70	!!
	71	!!--------------------------------------------------------------------
	72	INTEGER :: jh, nhan, jk, ji
[4147]	73	INTEGER :: ios ! Local integer output status for namelist read
[2956]	74
[3294]	75	NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
[2956]	76	!!----------------------------------------------------------------------
	77
[3294]	78	IF(lwp) THEN
	79	WRITE(numout,*)
	80	WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
	81	WRITE(numout,*) '~~~~~~~ '
	82	ENDIF
	83	!
[7646]	84	IF( .NOT. ln_tide ) CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
	85	!
[3294]	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
[6140]	137	IF (TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN
[2956]	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
[6140]	196	ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztempsshn(ji,jj)ssmask (ji,jj)
	197	ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztempun_b(ji,jj)ssumask(ji,jj)
	198	ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztempvn_b(ji,jj)ssvmask(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
[7910]	227	REAL(wp), DIMENSION(jpi,jpj,jpmax_harmo,2) :: ana_amp
[3294]	228	!!--------------------------------------------------------------------
[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)
[6140]	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)
[6140]	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	!
	367	END SUBROUTINE dia_harm_end
[2956]	368
[4292]	369
[3294]	370	SUBROUTINE dia_wri_harm
[2956]	371	!!--------------------------------------------------------------------
	372	!! * ROUTINE dia_wri_harm *
	373	!!
	374	!! ** Purpose : Write tidal harmonic analysis results in a netcdf file
	375	!!--------------------------------------------------------------------
	376	CHARACTER(LEN=lc) :: cltext
	377	CHARACTER(LEN=lc) :: &
[3294]	378	cdfile_name_T , & ! name of the file created (T-points)
	379	cdfile_name_U , & ! name of the file created (U-points)
	380	cdfile_name_V ! name of the file created (V-points)
	381	INTEGER :: jh
[2956]	382	!!----------------------------------------------------------------------
	383
	384	IF(lwp) WRITE(numout,*) ' '
[3294]	385	IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
[2956]	386	IF(lwp) WRITE(numout,*) ' '
	387
	388	! A) Elevation
	389	!/////////////
	390	!
	391	DO jh = 1, nb_ana
	392	CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
	393	CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
	394	END DO
	395
	396	! B) ubar
	397	!/////////
	398	!
	399	DO jh = 1, nb_ana
	400	CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
	401	CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) )
	402	END DO
	403
	404	! C) vbar
	405	!/////////
	406	!
	407	DO jh = 1, nb_ana
[4683]	408	CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh ) )
	409	CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
[2956]	410	END DO
[4683]	411	!
[3294]	412	END SUBROUTINE dia_wri_harm
[2956]	413
[4292]	414
[3294]	415	SUBROUTINE SUR_DETERMINE(init)
[4683]	416	!!---------------------------------------------------------------------------------
	417	!! * ROUTINE SUR_DETERMINE *
	418	!!
	419	!!
	420	!!
	421	!!---------------------------------------------------------------------------------
	422	INTEGER, INTENT(in) :: init
	423	!
	424	INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
	425	REAL(wp) :: zval1, zval2, zx1
[7910]	426	REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2
	427	INTEGER , DIMENSION(jpincomax) :: ipos2, ipivot
[4683]	428	!---------------------------------------------------------------------------------
[3294]	429
[4683]	430	IF( init == 1 ) THEN
	431	IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
	432	IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
	433	!
	434	ztmp3(:,:) = 0._wp
	435	!
	436	DO jk1_sd = 1, nsparse
	437	DO jk2_sd = 1, nsparse
	438	nisparse(jk2_sd) = nisparse(jk2_sd)
	439	njsparse(jk2_sd) = njsparse(jk2_sd)
	440	IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
	441	ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) &
	442	& + valuesparse(jk1_sd)*valuesparse(jk2_sd)
	443	ENDIF
	444	END DO
[4292]	445	END DO
[4683]	446	!
	447	DO jj_sd = 1 ,ninco
	448	ipos1(jj_sd) = jj_sd
	449	ipos2(jj_sd) = jj_sd
	450	END DO
	451	!
	452	DO ji_sd = 1 , ninco
	453	!
	454	!find greatest non-zero pivot:
	455	zval1 = ABS(ztmp3(ji_sd,ji_sd))
	456	!
	457	ipivot(ji_sd) = ji_sd
	458	DO jj_sd = ji_sd, ninco
	459	zval2 = ABS(ztmp3(ji_sd,jj_sd))
[6140]	460	IF( zval2 >= zval1 )THEN
[4683]	461	ipivot(ji_sd) = jj_sd
	462	zval1 = zval2
	463	ENDIF
	464	END DO
	465	!
	466	DO ji1_sd = 1, ninco
	467	zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd)
	468	zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd))
	469	ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd)
	470	ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
	471	END DO
	472	!
	473	ipos2(ji_sd) = ipos1(ipivot(ji_sd))
	474	ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
	475	ipos1(ji_sd) = ipos2(ji_sd)
	476	ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
	477	zpivot(ji_sd) = ztmp3(ji_sd,ji_sd)
	478	DO jj_sd = 1, ninco
	479	ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
	480	END DO
	481	!
	482	DO ji2_sd = ji_sd+1, ninco
	483	zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
	484	DO jj_sd=1,ninco
	485	ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
	486	END DO
	487	END DO
	488	!
	489	END DO
	490	!
	491	ENDIF ! End init==1
[3294]	492
[4683]	493	DO ji_sd = 1, ninco
	494	ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
[3294]	495	DO ji2_sd = ji_sd+1, ninco
[4683]	496	ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
	497	END DO
	498	END DO
[3294]	499
[4683]	500	!system solving:
	501	ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
	502	ji_sd = ninco
	503	DO ji_sd = ninco-1, 1, -1
	504	zx1 = 0._wp
	505	DO jj_sd = ji_sd+1, ninco
	506	zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
	507	END DO
	508	ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
	509	END DO
[3294]	510
[4683]	511	DO jj_sd =1, ninco
	512	ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
	513	END DO
[3294]	514
[4683]	515	!
	516	END SUBROUTINE SUR_DETERMINE
[3294]	517
[2956]	518	#else
	519	!!----------------------------------------------------------------------
	520	!! Default case : Empty module
	521	!!----------------------------------------------------------------------
	522	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .FALSE.
	523	CONTAINS
	524	SUBROUTINE dia_harm ( kt ) ! Empty routine
	525	INTEGER, INTENT( IN ) :: kt
	526	WRITE(,) 'dia_harm: you should not have seen this print'
	527	END SUBROUTINE dia_harm
[4292]	528	#endif
[2956]	529
	530	!!======================================================================
	531	END MODULE diaharm

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 7910

Download in other formats: