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

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

Last change on this file since 12685 was 12076, checked in by smueller, 5 years ago

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