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 branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 4634

Last change on this file since 4634 was 4634, checked in by clem, 10 years ago

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