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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago

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