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