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 @ 3038

Last change on this file since 3038 was 3038, checked in by cbricaud, 12 years ago

Add comment line
File size: 15.1 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  :   Initialize frequency array and  nodal factor for nit000_han
62      !!
63      !! History :
64      !!   9.0  O. Le Galloudec and J. Chanut (Original)
65      !!--------------------------------------------------------------------
66      !! * Local declarations
67      INTEGER :: jh, nhan, jk, ji
68      NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, nb_ana, tname
69
70      !!----------------------------------------------------------------------
71
72      ! Read namelist parameters:
73      ! -------------------------
74      REWIND ( numnam )
75      READ   ( numnam, nam_diaharm )
76
77      IF(lwp) WRITE(numout,*)
78      IF(lwp) WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
79      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
80
81      IF(lwp) WRITE(numout,*) 'First time step used for analysis:  nit000_han= ', nit000_han
82      IF(lwp) WRITE(numout,*) 'Last  time step used for analysis:  nitend_han= ', nitend_han
83      IF(lwp) WRITE(numout,*) 'Time step frequency for harmonic analysis:  nstep_han= ', nstep_han
84
85      IF (nb_ana > nb_harmo_max) THEN
86        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
87                                & nb_ana must be lower than nb_harmo_max, stop'
88        IF(lwp) WRITE(numout,*) 'nb_harmo_max= ', nb_harmo_max
89        nstop = nstop + 1
90      ENDIF
91
92      ! Basic checks on harmonic analysis time window:
93      ! ----------------------------------------------
94      IF (nit000 > nit000_han) THEN
95        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
96                                & nit000_han must be greater than nit000, stop'
97        IF(lwp) WRITE(numout,*) 'restart capability not implemented'
98        nstop = nstop + 1
99      ENDIF
100      IF (nitend < nitend_han) THEN
101        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
102                                & nitend_han must be lower than nitend, stop'
103        IF(lwp) WRITE(numout,*) 'restart capability not implemented'
104        nstop = nstop + 1
105      ENDIF
106
107      IF (MOD(nitend_han-nit000_han+1,nstep_han).NE.0) THEN
108        IF(lwp) WRITE(numout,*) ' E R R O R : dia_harm_init. &
109                                & analysis time span must be a multiple of nstep_han, stop'
110        nstop = nstop + 1
111      END IF
112
113      CALL tide_init_Wave
114
115      ALLOCATE(name    (nb_ana))
116      DO jk=1,nb_ana
117       DO ji=1,jpmax_harmo
118          IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN
119             name(jk) = ji
120             EXIT
121          END IF
122       END DO
123      END DO
124
125      ! Initialize frequency array:
126      ! ---------------------------
127      ALLOCATE(ana_freq(nb_ana))
128      ALLOCATE(vt      (nb_ana))
129      ALLOCATE(ut      (nb_ana))
130      ALLOCATE(ft      (nb_ana))
131
132      CALL tide_harmo(ana_freq, vt, ut , ft, name ,nb_ana)
133
134      IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
135
136      DO jh = 1, nb_ana
137        IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
138      END DO
139
140      ! Initialize temporary arrays:
141      ! ----------------------------
142      ALLOCATE( ana_temp(jpi,jpj,nb_ana*2,3))
143      ana_temp(:,:,:,:) = 0.e0
144
145   END SUBROUTINE dia_harm_init
146   
147   SUBROUTINE dia_harm ( kt )
148      !!----------------------------------------------------------------------
149      !!                 ***  ROUTINE dia_harm  ***
150      !!----------------------------------------------------------------------
151      !!         
152      !! ** Purpose :   Tidal harmonic analysis main routine
153      !!
154      !! ** Action  :   Sums ssh/u/v over time analysis [nit000_han,nitend_han]
155      !!
156      !! History :
157      !!   9.0  O. Le Galloudec and J. Chanut (Original)
158      !!--------------------------------------------------------------------
159      !! * Argument:
160      INTEGER, INTENT( IN ) :: kt
161
162      !! * Local declarations
163      INTEGER :: ji, jj, jh, jc, nhc
164      REAL(wp) :: ztime, ztemp
165
166      IF ( kt .EQ. nit000 ) CALL dia_harm_init
167
168      IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. &
169           (MOD(kt,nstep_han).EQ.0) ) THEN
170
171        ztime = kt*rdt 
172       
173        nhc = 0
174        DO jh = 1,nb_ana
175          DO jc = 1,2
176            nhc = nhc+1
177            ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
178                    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
179
180            DO jj = 1,jpj
181              DO ji = 1,jpi
182                ! Elevation
183                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1)                &
184                                        + ztemp*sshn(ji,jj)*tmask(ji,jj,1)       
185#if defined key_dynspg_ts
186                ! ubar
187                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2)                &
188                                        + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)
189                ! vbar
190                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3)                &
191                                        + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)
192#endif
193              END DO
194            END DO
195
196          END DO
197        END DO
198       
199      END IF
200
201      IF ( kt .EQ. nitend_han ) CALL dia_harm_end
202
203 
204   END SUBROUTINE dia_harm
205
206   SUBROUTINE dia_harm_end
207      !!----------------------------------------------------------------------
208      !!                 ***  ROUTINE diaharm_end  ***
209      !!----------------------------------------------------------------------
210      !!         
211      !! ** Purpose :  Compute the Real and Imaginary part of tidal constituents
212      !!
213      !! ** Action  :  Decompose the signal on the harmonic constituents
214      !!
215      !! History :
216      !!   9.0  O. Le Galloudec and J. Chanut (Original)
217      !!--------------------------------------------------------------------
218
219      !! * Local declarations
220      INTEGER :: ji, jj, jh, jc, jn, nhan, jl
221      INTEGER :: ksp, kun, keq
222      REAL(wp) :: ztime, ztime_ini, ztime_end
223      REAL(wp) :: X1,X2
224      REAL(wp), DIMENSION(jpi,jpj,nb_harmo_max,2) :: ana_amp
225
226
227        IF(lwp) WRITE(numout,*)
228        IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
229        IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
230
231        ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis
232        ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis
233        nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
234
235        NBINCO = 2*nb_ana
236
237        ksp = 0
238        keq = 0       
239        DO jn = 1, nhan
240          ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1)
241          keq = keq + 1
242          kun = 0
243          DO jh = 1,nb_ana
244            DO jc = 1,2
245              kun = kun + 1
246              ksp = ksp + 1
247              ISPARSE(ksp) = keq
248              JSPARSE(ksp) = kun
249              SPARSEVALUE(ksp)= &
250                 +(     MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
251                   +(1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
252            END DO
253          END DO
254        END DO
255
256        NBSPARSE=ksp
257
258        ! Elevation:
259        DO jj = 1, jpj
260          DO ji = 1, jpi
261            ! Fill input array
262            kun=0
263            DO jh = 1,nb_ana
264              DO jc = 1,2
265                kun = kun + 1
266                TAB4(kun)=ana_temp(ji,jj,kun,1)
267              ENDDO
268            ENDDO
269
270            CALL SUR_DETERMINE(jj)
271
272            ! Fill output array
273            DO jh = 1, nb_ana
274              ana_amp(ji,jj,jh,1)=TAB7((jh-1)*2+1)
275              ana_amp(ji,jj,jh,2)=TAB7((jh-1)*2+2)
276            END DO
277          END DO
278        END DO
279
280        ALLOCATE(out_eta(jpi,jpj,2*nb_ana))
281        ALLOCATE(out_u  (jpi,jpj,2*nb_ana))
282        ALLOCATE(out_v  (jpi,jpj,2*nb_ana))
283
284
285        DO jj = 1, jpj
286          DO ji = 1, jpi
287            DO jh = 1, nb_ana 
288                X1=ana_amp(ji,jj,jh,1)
289                X2=-ana_amp(ji,jj,jh,2)
290                out_eta(ji,jj,jh)=X1 * tmask(ji,jj,1)
291                out_eta(ji,jj,nb_ana+jh)=X2 * tmask(ji,jj,1)
292            ENDDO
293          ENDDO
294        ENDDO
295
296        ! ubar:
297        DO jj = 1, jpj
298          DO ji = 1, jpi
299            ! Fill input array
300            kun=0
301            DO jh = 1,nb_ana
302              DO jc = 1,2
303                kun = kun + 1
304                TAB4(kun)=ana_temp(ji,jj,kun,2)
305              ENDDO
306            ENDDO
307
308            CALL SUR_DETERMINE(jj+1)
309
310            ! Fill output array
311            DO jh = 1, nb_ana
312              ana_amp(ji,jj,jh,1)=TAB7((jh-1)*2+1)
313              ana_amp(ji,jj,jh,2)=TAB7((jh-1)*2+2)
314            END DO
315
316          END DO
317        END DO
318
319        DO jj = 1, jpj
320          DO ji = 1, jpi
321            DO jh = 1, nb_ana 
322                X1=ana_amp(ji,jj,jh,1)
323                X2=-ana_amp(ji,jj,jh,2)
324                 out_u(ji,jj,jh) = X1 * umask(ji,jj,1)
325                 out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1)
326            ENDDO
327          ENDDO
328        ENDDO
329
330        ! vbar:
331        DO jj = 1, jpj
332          DO ji = 1, jpi
333              ! Fill input array
334              kun=0
335              DO jh = 1,nb_ana
336                DO jc = 1,2
337                  kun = kun + 1
338                  TAB4(kun)=ana_temp(ji,jj,kun,3)
339                ENDDO
340              ENDDO
341
342              CALL SUR_DETERMINE(jj+1)
343
344              ! Fill output array
345              DO jh = 1, nb_ana
346                ana_amp(ji,jj,jh,1)=TAB7((jh-1)*2+1)
347                ana_amp(ji,jj,jh,2)=TAB7((jh-1)*2+2)
348              END DO
349
350          END DO
351        END DO
352
353        DO jj = 1, jpj
354          DO ji = 1, jpi
355            DO jh = 1, nb_ana 
356                X1=ana_amp(ji,jj,jh,1)
357                X2=-ana_amp(ji,jj,jh,2)
358                 out_v(ji,jj,jh)=X1 * vmask(ji,jj,1)
359                 out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1)
360            ENDDO
361          ENDDO
362        ENDDO
363
364        CALL dia_wri_harm ! Write results in files
365
366  END SUBROUTINE dia_harm_end
367
368  SUBROUTINE dia_wri_harm
369      !!--------------------------------------------------------------------
370      !!                 ***  ROUTINE dia_wri_harm  ***
371      !!--------------------------------------------------------------------
372      !!         
373      !! ** Purpose : Write tidal harmonic analysis results in a netcdf file
374      !!
375      !!
376      !! History :
377      !!   9.0  O. Le Galloudec and J. Chanut (Original)
378      !!--------------------------------------------------------------------
379
380      !! * Local declarations
381      CHARACTER(LEN=lc) :: cltext
382      CHARACTER(LEN=lc) ::   &
383         cdfile_name_T    ,   & ! name of the file created (T-points)
384         cdfile_name_U    ,   & ! name of the file created (U-points)
385         cdfile_name_V          ! name of the file created (V-points)
386      INTEGER  ::   jh       
387      !!----------------------------------------------------------------------
388
389#if defined key_dimgout
390      cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'
391      cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'
392      cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'
393#endif
394
395      IF(lwp) WRITE(numout,*) '  '
396      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 
397#if defined key_dimgout
398      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T)
399      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_U)
400      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_V)
401#endif
402      IF(lwp) WRITE(numout,*) '  '
403
404      ! A) Elevation
405      !/////////////
406      !
407#if defined key_dimgout
408      cltext='Elevation amplitude and phase'
409      CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')
410#else
411      DO jh = 1, nb_ana
412      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
413      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
414      END DO
415#endif
416
417      ! B) ubar
418      !/////////
419      !
420#if defined key_dimgout
421      cltext='ubar amplitude and phase'
422      CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')
423#else
424      DO jh = 1, nb_ana
425      CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
426      CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) )
427      END DO
428#endif
429
430      ! C) vbar
431      !/////////
432      !
433#if defined key_dimgout
434      cltext='vbar amplitude and phase'
435      CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')
436#else
437      DO jh = 1, nb_ana
438      CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh) )
439      CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,nb_ana+jh) )
440      END DO
441#endif
442
443  END SUBROUTINE dia_wri_harm
444
445#else
446   !!----------------------------------------------------------------------
447   !!   Default case :   Empty module
448   !!----------------------------------------------------------------------
449   LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE.
450CONTAINS
451
452   SUBROUTINE dia_harm ( kt )     ! Empty routine
453      INTEGER, INTENT( IN ) :: kt 
454      WRITE(*,*) 'dia_harm: you should not have seen this print'
455   END SUBROUTINE dia_harm
456
457
458#endif
459   !!======================================================================
460END MODULE diaharm
Note: See TracBrowser for help on using the repository browser.