Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

diaharm.F90

Last change on this file was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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

Download in other formats: