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

source: branches/2011/dev_r2787_MERCATOR2_tidalharm/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 2964

Last change on this file since 2964 was 2964, checked in by cbricaud, 13 years ago
File size: 14.9 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 surdetermine
18   USE daymod
19   USE tide_mod
20   USE iom 
21
22   IMPLICIT NONE
23   PRIVATE
24   
25   INTEGER, PARAMETER :: nb_harmo_max=9
26
27   LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm  = .TRUE.
28
29   INTEGER ::      & !! namelist variables
30       nit000_han=1, & ! First time step used for harmonic analysis
31       nitend_han=1, & ! Last time step used for harmonic analysis
32       nstep_han=1,  & ! Time step frequency for harmonic analysis
33       nb_ana          ! Number of harmonics to analyse
34
35   CHARACTER (LEN=4), DIMENSION(nb_harmo_max) ::   &
36       tname         ! Names of tidal constituents ('M2', 'K1',...)
37
38   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
39   REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, vt, ut, ft
40   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta, &
41              out_u, &
42                 out_v
43   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: name
44
45!! * Routine accessibility
46   PUBLIC  dia_harm    ! routine called by step.F90
47
48   !!---------------------------------------------------------------------------------
49   !! 
50   !!---------------------------------------------------------------------------------
51
52CONTAINS
53
54   SUBROUTINE dia_harm_init 
55      !!----------------------------------------------------------------------
56      !!                 ***  ROUTINE dia_harm_init  ***
57      !!----------------------------------------------------------------------
58      !!         
59      !! ** Purpose :   Initialization of tidal harmonic analysis
60      !!
61      !! ** Method  :   
62      !!
63      !! ** Input   : 
64      !!
65      !! History :
66      !!   9.0  O. Le Galloudec and J. Chanut (Original)
67      !!--------------------------------------------------------------------
68      !! * Local declarations
69      INTEGER :: jh, nhan, jk, ji
70      NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, nb_ana, tname
71
72      !!----------------------------------------------------------------------
73
74      ! Read namelist parameters:
75      ! -------------------------
76      REWIND ( numnam )
77      READ   ( numnam, nam_diaharm )
78
79      IF(lwp) WRITE(numout,*)
80      IF(lwp) WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
81      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
82
83      IF(lwp) WRITE(numout,*) 'First time step used for analysis:  nit000_han= ', nit000_han
84      IF(lwp) WRITE(numout,*) 'Last  time step used for analysis:  nitend_han= ', nitend_han
85      IF(lwp) WRITE(numout,*) 'Time step frequency for harmonic analysis:  nstep_han= ', nstep_han
86
87      IF (nb_ana > nb_harmo_max) THEN
88        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
89                                & nb_ana must be lower than nb_harmo_max, stop'
90        IF(lwp) WRITE(numout,*) 'nb_harmo_max= ', nb_harmo_max
91        nstop = nstop + 1
92      ENDIF
93
94      ! Basic checks on harmonic analysis time window:
95      ! ----------------------------------------------
96      IF (nit000 > nit000_han) THEN
97        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
98                                & nit000_han must be greater than nit000, stop'
99        IF(lwp) WRITE(numout,*) 'restart capability not implemented'
100        nstop = nstop + 1
101      ENDIF
102      IF (nitend < nitend_han) THEN
103        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
104                                & nitend_han must be lower than nitend, stop'
105        IF(lwp) WRITE(numout,*) 'restart capability not implemented'
106        nstop = nstop + 1
107      ENDIF
108
109      IF (MOD(nitend_han-nit000_han+1,nstep_han).NE.0) THEN
110        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
111                                & analysis time span must be a multiple of nstep_han, stop'
112        nstop = nstop + 1
113      END IF
114
115      CALL tide_init_Wave
116
117      ALLOCATE(name    (nb_ana))
118      DO jk=1,nb_ana
119       DO ji=1,jpmax_harmo
120          IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN
121             name(jk) = ji
122             EXIT
123          END IF
124       END DO
125      END DO
126
127      ! Initialize frequency array:
128      ! ---------------------------
129      ALLOCATE(ana_freq(nb_ana))
130      ALLOCATE(vt      (nb_ana))
131      ALLOCATE(ut      (nb_ana))
132      ALLOCATE(ft      (nb_ana))
133
134      CALL tide_harmo(ana_freq, vt, ut , ft, name ,nb_ana)
135
136      IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
137
138      DO jh = 1, nb_ana
139        IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
140      END DO
141
142      ! Initialize temporary arrays:
143      ! ----------------------------
144      ALLOCATE( ana_temp(jpi,jpj,nb_ana*2,3))
145      ana_temp(:,:,:,:) = 0.e0
146
147   END SUBROUTINE dia_harm_init
148   
149   SUBROUTINE dia_harm ( kt )
150      !!----------------------------------------------------------------------
151      !!                 ***  ROUTINE dia_harm  ***
152      !!----------------------------------------------------------------------
153      !!         
154      !! ** Purpose :   Tidal harmonic analysis main routine
155      !!
156      !! ** Method  :   
157      !!
158      !! ** Input   : 
159      !!
160      !! History :
161      !!   9.0  O. Le Galloudec and J. Chanut (Original)
162      !!--------------------------------------------------------------------
163      !! * Argument:
164      INTEGER, INTENT( IN ) :: kt
165
166      !! * Local declarations
167      INTEGER :: ji, jj, jh, jc, nhc
168      REAL(wp) :: ztime, ztemp
169
170      IF ( kt .EQ. nit000 ) CALL dia_harm_init
171
172      IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. &
173           (MOD(kt,nstep_han).EQ.0) ) THEN
174
175        ztime = kt*rdt 
176       
177        nhc = 0
178        DO jh = 1,nb_ana
179          DO jc = 1,2
180            nhc = nhc+1
181            ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
182                    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
183
184            DO jj = 1,jpj
185              DO ji = 1,jpi
186                ! Elevation
187                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1)                &
188                                        + ztemp*sshn(ji,jj)*tmask(ji,jj,1)       
189#if defined key_dynspg_ts
190                ! ubar
191                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2)                &
192                                        + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)
193                ! vbar
194                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3)                &
195                                        + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)
196#endif
197              END DO
198            END DO
199
200          END DO
201        END DO
202       
203      END IF
204
205      IF ( kt .EQ. nitend_han ) CALL dia_harm_end
206
207 
208   END SUBROUTINE dia_harm
209
210   SUBROUTINE dia_harm_end
211      !!----------------------------------------------------------------------
212      !!                 ***  ROUTINE diaharm_end  ***
213      !!----------------------------------------------------------------------
214      !!         
215      !! ** Purpose :
216      !!
217      !! ** Method  :   
218      !!
219      !! ** Input   : 
220      !!
221      !! History :
222      !!   9.0  O. Le Galloudec and J. Chanut (Original)
223      !!--------------------------------------------------------------------
224
225      !! * Local declarations
226      INTEGER :: ji, jj, jh, jc, jn, nhan, jl
227      INTEGER :: ksp, kun, keq
228      REAL(wp) :: ztime, ztime_ini, ztime_end
229      REAL(wp) :: X1,X2
230      REAL(wp), DIMENSION(jpi,jpj,nb_harmo_max,2) :: ana_amp
231
232
233        IF(lwp) WRITE(numout,*)
234        IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
235        IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
236
237        ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis
238        ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis
239        nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
240
241        NBINCO = 2*nb_ana
242
243        ksp = 0
244        keq = 0       
245        DO jn = 1, nhan
246          ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1)
247          keq = keq + 1
248          kun = 0
249          DO jh = 1,nb_ana
250            DO jc = 1,2
251              kun = kun + 1
252              ksp = ksp + 1
253              ISPARSE(ksp) = keq
254              JSPARSE(ksp) = kun
255              SPARSEVALUE(ksp)= &
256                 +(     MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
257                   +(1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
258            END DO
259          END DO
260        END DO
261
262        NBSPARSE=ksp
263
264        ! Elevation:
265        DO jj = 1, jpj
266          DO ji = 1, jpi
267            ! Fill input array
268            kun=0
269            DO jh = 1,nb_ana
270              DO jc = 1,2
271                kun = kun + 1
272                TAB4(kun)=ana_temp(ji,jj,kun,1)
273              ENDDO
274            ENDDO
275
276            CALL SUR_DETERMINE(jj)
277
278            ! Fill output array
279            DO jh = 1, nb_ana
280              ana_amp(ji,jj,jh,1)=TAB7((jh-1)*2+1)
281              ana_amp(ji,jj,jh,2)=TAB7((jh-1)*2+2)
282            END DO
283          END DO
284        END DO
285
286        ALLOCATE(out_eta(jpi,jpj,2*nb_ana))
287        ALLOCATE(out_u  (jpi,jpj,2*nb_ana))
288        ALLOCATE(out_v  (jpi,jpj,2*nb_ana))
289
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,nb_ana+jh)=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                TAB4(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)=TAB7((jh-1)*2+1)
319              ana_amp(ji,jj,jh,2)=TAB7((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                  TAB4(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)=TAB7((jh-1)*2+1)
353                ana_amp(ji,jj,jh,2)=TAB7((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
372  END SUBROUTINE dia_harm_end
373
374  SUBROUTINE dia_wri_harm
375      !!--------------------------------------------------------------------
376      !!                 ***  ROUTINE dia_wri_harm  ***
377      !!--------------------------------------------------------------------
378      !!         
379      !! ** Purpose : Write tidal harmonic analysis results in a netcdf file
380      !!
381      !!
382      !! History :
383      !!   9.0  O. Le Galloudec and J. Chanut (Original)
384      !!--------------------------------------------------------------------
385
386      !! * Local declarations
387      CHARACTER(LEN=lc) :: cltext
388      CHARACTER(LEN=lc) ::   &
389         cdfile_name_T    ,   & ! name of the file created (T-points)
390         cdfile_name_U    ,   & ! name of the file created (U-points)
391         cdfile_name_V          ! name of the file created (V-points)
392      INTEGER  ::   jh       
393      !!----------------------------------------------------------------------
394
395#if defined key_dimgout
396      cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'
397      cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'
398      cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'
399#endif
400
401      IF(lwp) WRITE(numout,*) '  '
402      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 
403#if defined key_dimgout
404      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T)
405      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_U)
406      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_V)
407#endif
408      IF(lwp) WRITE(numout,*) '  '
409
410      ! A) Elevation
411      !/////////////
412      !
413#if defined key_dimgout
414      cltext='Elevation amplitude and phase'
415      CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')
416#else
417      DO jh = 1, nb_ana
418      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
419      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
420      END DO
421#endif
422
423      ! B) ubar
424      !/////////
425      !
426#if defined key_dimgout
427      cltext='ubar amplitude and phase'
428      CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')
429#else
430      DO jh = 1, nb_ana
431      CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
432      CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) )
433      END DO
434#endif
435
436      ! C) vbar
437      !/////////
438      !
439#if defined key_dimgout
440      cltext='vbar amplitude and phase'
441      CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')
442#else
443      DO jh = 1, nb_ana
444      CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh) )
445      CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,nb_ana+jh) )
446      END DO
447#endif
448
449  END SUBROUTINE dia_wri_harm
450
451#else
452   !!----------------------------------------------------------------------
453   !!   Default case :   Empty module
454   !!----------------------------------------------------------------------
455   LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE.
456CONTAINS
457
458   SUBROUTINE dia_harm ( kt )     ! Empty routine
459      INTEGER, INTENT( IN ) :: kt 
460      WRITE(*,*) 'dia_harm: you should not have seen this print'
461   END SUBROUTINE dia_harm
462
463
464#endif
465   !!======================================================================
466END MODULE diaharm
Note: See TracBrowser for help on using the repository browser.