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 NEMO/branches/UKMO/NEMO_4.0.4_mirror/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_mirror/src/OCE/DIA/diaharm.F90 @ 14075

Last change on this file since 14075 was 14075, checked in by cguiavarch, 22 months ago

UKMO/NEMO_4.0.4_mirror : Remove SVN keywords.

File size: 17.4 KB
Line 
1MODULE diaharm 
2   !!======================================================================
3   !!                       ***  MODULE  diaharm  ***
4   !! Harmonic analysis of tidal constituents
5   !!======================================================================
6   !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code
7   !!----------------------------------------------------------------------
8   USE oce             ! ocean dynamics and tracers variables
9   USE dom_oce         ! ocean space and time domain
10   USE phycst
11   USE daymod
12   USE tide_mod
13   USE sbctide         ! Tidal forcing or not
14   !
15   USE in_out_manager  ! I/O units
16   USE iom             ! I/0 library
17   USE ioipsl          ! NetCDF IPSL library
18   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
19   USE timing          ! preformance summary
20   USE lib_mpp           ! MPP library
21
22   IMPLICIT NONE
23   PRIVATE
24   
25   INTEGER, PARAMETER :: jpincomax    = 2.*jpmax_harmo
26   INTEGER, PARAMETER :: jpdimsparse  = jpincomax*366*24*2 ! 30min for a 1yr-long run
27
28   !                         !!** namelist variables **
29   LOGICAL, PUBLIC ::   ln_diaharm    ! Choose tidal harmonic output or not
30   INTEGER         ::   nit000_han    ! First time step used for harmonic analysis
31   INTEGER         ::   nitend_han    ! Last time step used for harmonic analysis
32   INTEGER         ::   nstep_han     ! Time step frequency for harmonic analysis
33   INTEGER         ::   nb_ana        ! Number of harmonics to analyse
34
35   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name
36   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ana_temp
37   REAL(wp), ALLOCATABLE, DIMENSION(:)       ::   ana_freq, ut, vt, ft
38   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   out_eta, out_u, out_v
39
40   INTEGER  ::   ninco, nsparse
41   REAL(wp) ::   z1_tmp3
42   INTEGER ,       DIMENSION(jpdimsparse)         ::   njsparse, nisparse
43   INTEGER , SAVE, DIMENSION(jpincomax)           ::   ipos1
44   REAL(wp),       DIMENSION(jpdimsparse)         ::   valuesparse
45   REAL(wp),       DIMENSION(jpincomax)           ::   ztmp4 , ztmp7, z1_pivot
46   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) ::   ztmp3 , zpilier
47
48   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
49
50   PUBLIC   dia_harm        ! routine called by step.F90
51   PUBLIC   dia_harm_init   ! routine called by nemogcm.F90
52
53   !!----------------------------------------------------------------------
54   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
55   !! $Id$
56   !! Software governed by the CeCILL license (see ./LICENSE)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE dia_harm_init 
61      !!----------------------------------------------------------------------
62      !!                 ***  ROUTINE dia_harm_init  ***
63      !!         
64      !! ** Purpose :   Initialization of tidal harmonic analysis
65      !!
66      !! ** Method  :   Initialize frequency array and  nodal factor for nit000_han
67      !!
68      !!--------------------------------------------------------------------
69      INTEGER ::   jh, nhan, ji
70      INTEGER ::   ios                 ! Local integer output status for namelist read
71
72      NAMELIST/nam_diaharm/ ln_diaharm, nit000_han, nitend_han, nstep_han, tname
73      !!----------------------------------------------------------------------
74
75      IF(lwp) THEN
76         WRITE(numout,*)
77         WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
78         WRITE(numout,*) '~~~~~~~~~~~~~ '
79      ENDIF
80      !
81      REWIND( numnam_ref )              ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
82      READ  ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
83901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_diaharm in reference namelist' )
84      REWIND( numnam_cfg )              ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis
85      READ  ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 )
86902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist' )
87      IF(lwm) WRITE ( numond, nam_diaharm )
88      !
89      IF(lwp) THEN
90         WRITE(numout,*) 'Tidal diagnostics = ', ln_diaharm
91         WRITE(numout,*) '   First time step used for analysis:         nit000_han= ', nit000_han
92         WRITE(numout,*) '   Last  time step used for analysis:         nitend_han= ', nitend_han
93         WRITE(numout,*) '   Time step frequency for harmonic analysis: nstep_han = ', nstep_han
94      ENDIF
95
96      IF( ln_diaharm .AND. .NOT.ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
97
98      IF( ln_diaharm ) THEN
99
100         CALL tide_init_Wave
101         !
102         ! Basic checks on harmonic analysis time window:
103         ! ----------------------------------------------
104         IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   &
105            &                                       ' restart capability not implemented' )
106         IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   &
107            &                                       'restart capability not implemented' )
108
109         IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   &
110            &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
111         !
112         nb_ana = 0
113         DO jh=1,jpmax_harmo
114            DO ji=1,jpmax_harmo
115               IF(TRIM(tname(jh)) == Wave(ji)%cname_tide) THEN
116                  nb_ana=nb_ana+1
117               ENDIF
118            END DO
119         END DO
120         !
121         IF(lwp) THEN
122            WRITE(numout,*) '        Namelist nam_diaharm'
123            WRITE(numout,*) '        nb_ana    = ', nb_ana
124            CALL flush(numout)
125         ENDIF
126         !
127         IF (nb_ana > jpmax_harmo) THEN
128            WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo'
129            WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo
130            CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 )
131         ENDIF
132
133         ALLOCATE(name(nb_ana))
134         DO jh=1,nb_ana
135            DO ji=1,jpmax_harmo
136               IF (TRIM(tname(jh)) ==  Wave(ji)%cname_tide) THEN
137                  name(jh) = ji
138                  EXIT
139               END IF
140            END DO
141         END DO
142
143         ! Initialize frequency array:
144         ! ---------------------------
145         ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
146
147         CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
148
149         IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
150
151         DO jh = 1, nb_ana
152            IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
153         END DO
154
155         ! Initialize temporary arrays:
156         ! ----------------------------
157         ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
158         ana_temp(:,:,:,:) = 0._wp
159
160      ENDIF
161
162   END SUBROUTINE dia_harm_init
163
164
165   SUBROUTINE dia_harm( kt )
166      !!----------------------------------------------------------------------
167      !!                 ***  ROUTINE dia_harm  ***
168      !!         
169      !! ** Purpose :   Tidal harmonic analysis main routine
170      !!
171      !! ** Action  :   Sums ssh/u/v over time analysis [nit000_han,nitend_han]
172      !!
173      !!--------------------------------------------------------------------
174      INTEGER, INTENT( in ) ::   kt
175      !
176      INTEGER  ::   ji, jj, jh, jc, nhc
177      REAL(wp) ::   ztime, ztemp
178      !!--------------------------------------------------------------------
179      IF( ln_timing )   CALL timing_start('dia_harm')
180      !
181      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
182         !
183         ztime = (kt-nit000+1) * rdt 
184         !
185         nhc = 0
186         DO jh = 1, nb_ana
187            DO jc = 1, 2
188               nhc = nhc+1
189               ztemp =  (   MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
190                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
191                  !
192               DO jj = 2, jpjm1
193                  DO ji = 2, jpim1
194                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp * sshn(ji,jj) * ssmask (ji,jj) ! elevation     
195                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp * un_b(ji,jj) * ssumask(ji,jj) ! u-vel
196                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp * vn_b(ji,jj) * ssvmask(ji,jj) ! v-vel
197                  END DO
198               END DO
199            END DO
200         END DO
201      END IF
202      !
203      IF( kt == nitend_han )   CALL dia_harm_end
204      !
205      IF( ln_timing )   CALL timing_stop('dia_harm')
206      !
207   END SUBROUTINE dia_harm
208
209
210   SUBROUTINE dia_harm_end
211      !!----------------------------------------------------------------------
212      !!                 ***  ROUTINE diaharm_end  ***
213      !!         
214      !! ** Purpose :  Compute the Real and Imaginary part of tidal constituents
215      !!
216      !! ** Action  :  Decompose the signal on the harmonic constituents
217      !!
218      !!--------------------------------------------------------------------
219      INTEGER  ::   ji, jj, jh, jc, jn, nhan
220      INTEGER  ::   ksp, kun, keq
221      REAL(wp) ::   ztime, ztime_ini, ztime_end, z1_han
222      !!--------------------------------------------------------------------
223      !
224      IF(lwp) WRITE(numout,*)
225      IF(lwp) WRITE(numout,*) 'dia_harm_end: kt=nitend_han: Perform harmonic analysis'
226      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
227     
228      ALLOCATE( out_eta(jpi,jpj,2*nb_ana), out_u(jpi,jpj,2*nb_ana), out_v(jpi,jpj,2*nb_ana) )
229
230      ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis
231      ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis
232      nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
233      z1_han = 1._wp / REAL(nhan-1) 
234     
235      ninco = 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 ) * z1_han
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               nisparse(ksp) = keq
248               njsparse(ksp) = kun
249               valuesparse(ksp) = (   MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))   &
250                  &             + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
251            END DO
252         END DO
253      END DO
254
255      nsparse = ksp
256
257      IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
258      IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
259
260      CALL SUR_DETERMINE_INIT
261
262      ! Elevation:
263      DO jj = 2, jpjm1
264         DO ji = 2, jpim1
265
266            ! Fill input array
267            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1)
268            CALL SUR_DETERMINE
269           
270            ! Fill output array
271            DO jh = 1, nb_ana
272               out_eta(ji,jj,jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj)
273               out_eta(ji,jj,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj)
274            END DO
275         END DO
276      END DO
277
278      ! ubar:
279      DO jj = 2, jpjm1
280         DO ji = 2, jpim1
281
282            ! Fill input array
283            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2)
284            CALL SUR_DETERMINE
285
286            ! Fill output array
287            DO jh = 1, nb_ana
288               out_u(ji,jj,       jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj)
289               out_u(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj)
290            END DO
291
292        END DO
293      END DO
294
295      ! vbar:
296      DO jj = 2, jpjm1
297         DO ji = 2, jpim1
298
299            ! Fill input array
300            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3)
301            CALL SUR_DETERMINE
302
303            ! Fill output array
304            DO jh = 1, nb_ana
305               out_v(ji,jj,       jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj)
306               out_v(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj)
307            END DO
308
309         END DO
310      END DO
311      !
312      ! clem: we could avoid this call if all the loops were from 1:jpi and 1:jpj
313      !       but I think this is the most efficient
314      CALL lbc_lnk_multi( 'dia_harm_end', out_eta, 'T', 1., out_u, 'U', -1. , out_v, 'V', -1. )
315      !
316      CALL dia_wri_harm ! Write results in files
317      !
318      DEALLOCATE( out_eta, out_u, out_v )
319      !
320   END SUBROUTINE dia_harm_end
321
322
323   SUBROUTINE dia_wri_harm
324      !!--------------------------------------------------------------------
325      !!                 ***  ROUTINE dia_wri_harm  ***
326      !!         
327      !! ** Purpose : Write tidal harmonic analysis results in a netcdf file
328      !!--------------------------------------------------------------------
329      INTEGER  ::   jh
330      !!----------------------------------------------------------------------
331
332      IF(lwp) WRITE(numout,*) '  '
333      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
334      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
335
336      ! A) Elevation
337      !/////////////
338      DO jh = 1, nb_ana
339      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
340      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,jh+nb_ana) )
341      END DO
342
343      ! B) ubar
344      !/////////
345      DO jh = 1, nb_ana
346      CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
347      CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,jh+nb_ana) )
348      END DO
349
350      ! C) vbar
351      !/////////
352      DO jh = 1, nb_ana
353         CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) )
354         CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
355      END DO
356      !
357   END SUBROUTINE dia_wri_harm
358
359
360   SUBROUTINE SUR_DETERMINE_INIT
361      !!---------------------------------------------------------------------------------
362      !!                      *** ROUTINE SUR_DETERMINE_INIT ***
363      !!       
364      !!---------------------------------------------------------------------------------
365      INTEGER                        :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd
366      INTEGER                        :: ipivot
367      REAL(wp)                       :: zval1, zval2, zcol1, zcol2
368      INTEGER , DIMENSION(jpincomax) :: ipos2
369      !!---------------------------------------------------------------------------------
370      !           
371      !
372      ztmp3(:,:) = 0._wp
373      !
374      DO jh1_sd = 1, nsparse
375         DO jh2_sd = 1, nsparse
376            IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN
377               ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  &
378                  &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd)
379            ENDIF
380         END DO
381      END DO
382      !
383      DO jj_sd = 1, ninco
384         ipos1(jj_sd) = jj_sd
385         ipos2(jj_sd) = jj_sd
386      END DO
387      !
388      DO ji_sd = 1, ninco
389         !
390         !find greatest non-zero pivot:
391         zval1 = ABS(ztmp3(ji_sd,ji_sd))
392         !
393         ipivot = ji_sd
394         DO jj_sd = ji_sd, ninco
395            zval2 = ABS(ztmp3(ji_sd,jj_sd))
396            IF( zval2 >= zval1 )THEN
397               ipivot = jj_sd
398               zval1  = zval2
399            ENDIF
400         END DO
401         !
402         DO ji1_sd = 1, ninco
403            zcol1                = ztmp3(ji1_sd,ji_sd)
404            zcol2                = ztmp3(ji1_sd,ipivot)
405            ztmp3(ji1_sd,ji_sd)  = zcol2
406            ztmp3(ji1_sd,ipivot) = zcol1
407         END DO
408         !
409         ipos2(ji_sd)    = ipos1(ipivot)
410         ipos2(ipivot)   = ipos1(ji_sd)
411         ipos1(ji_sd)    = ipos2(ji_sd)
412         ipos1(ipivot)   = ipos2(ipivot)
413         z1_pivot(ji_sd) = 1._wp / ztmp3(ji_sd,ji_sd)
414         DO jj_sd = 1, ninco
415            ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) * z1_pivot(ji_sd)
416         END DO
417         !
418         DO ji2_sd = ji_sd+1, ninco
419            zpilier(ji2_sd,ji_sd) = ztmp3(ji2_sd,ji_sd)
420            DO jj_sd=1,ninco
421               ztmp3(ji2_sd,jj_sd) = ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
422            END DO
423         END DO
424         !
425      END DO
426      !
427      z1_tmp3 = 1._wp / ztmp3(ninco,ninco)
428      !
429   END SUBROUTINE SUR_DETERMINE_INIT
430
431   
432   SUBROUTINE SUR_DETERMINE
433      !!---------------------------------------------------------------------------------
434      !!                      *** ROUTINE SUR_DETERMINE ***
435      !!   
436      !!---------------------------------------------------------------------------------
437      INTEGER                        :: ji_sd, jj_sd, ji1_sd, ji2_sd
438      REAL(wp)                       :: zx1
439      REAL(wp), DIMENSION(jpincomax) :: ztmpx
440      !!---------------------------------------------------------------------------------
441      !           
442      DO ji_sd = 1, ninco
443         ztmp4(ji_sd) = ztmp4(ji_sd) * z1_pivot(ji_sd)
444         DO ji2_sd = ji_sd+1, ninco
445            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
446         END DO
447      END DO
448
449      !system solving:
450      ztmpx(ninco) = ztmp4(ninco) * z1_tmp3
451      DO ji_sd = ninco-1, 1, -1
452         zx1 = 0._wp
453         DO jj_sd = ji_sd+1, ninco
454            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
455         END DO
456         ztmpx(ji_sd) = ztmp4(ji_sd) - zx1
457      END DO
458
459      DO jj_sd = 1, ninco
460         ztmp7(ipos1(jj_sd)) = ztmpx(jj_sd)
461      END DO
462      !
463   END SUBROUTINE SUR_DETERMINE
464
465   
466   !!======================================================================
467END MODULE diaharm
Note: See TracBrowser for help on using the repository browser.