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/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DIA/diaharm.F90 @ 12076

Last change on this file since 12076 was 12076, checked in by smueller, 5 years ago

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