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

source: branches/2013/dev_r3853_CNRS9_ConfSetting/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 3901

Last change on this file since 3901 was 3901, checked in by clevy, 11 years ago

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