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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 3194

Last change on this file since 3194 was 3186, checked in by smasson, 13 years ago

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