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

source: branches/2011/dev_LOCEAN_CMCC_INGV_MERCATOR_2011/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 3104

Last change on this file since 3104 was 3104, checked in by cetlod, 12 years ago

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