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_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/DIA/diaharm.F90 @ 11671

Last change on this file since 11671 was 11671, checked in by acc, 5 years ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Final, non-substantive changes to complete this branch. These changes remove all REWIND statements on the old namelist fortran units (now character variables for internal files). These changes have been left until last since they are easily repeated via a script and it may be preferable to use the previous revision for merge purposes and reapply these last changes separately. This branch has been fully SETTE tested.
  • Property svn:keywords set to Id
File size: 18.3 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*300*24
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   INTEGER ,       DIMENSION(jpdimsparse)         ::   njsparse, nisparse
42   INTEGER , SAVE, DIMENSION(jpincomax)           ::   ipos1
43   REAL(wp),       DIMENSION(jpdimsparse)         ::   valuesparse
44   REAL(wp),       DIMENSION(jpincomax)           ::   ztmp4 , ztmp7
45   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) ::   ztmp3 , zpilier
46   REAL(wp), SAVE, DIMENSION(jpincomax)           ::   zpivot
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      READ  ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
82901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_diaharm in reference namelist' )
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         CALL tide_init_Wave
99         !
100         ! Basic checks on harmonic analysis time window:
101         ! ----------------------------------------------
102         IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   &
103            &                                       ' restart capability not implemented' )
104         IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   &
105            &                                       'restart capability not implemented' )
106
107         IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   &
108            &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
109         !
110         nb_ana = 0
111         DO jh=1,jpmax_harmo
112            DO ji=1,jpmax_harmo
113               IF(TRIM(tname(jh)) == Wave(ji)%cname_tide) THEN
114                  nb_ana=nb_ana+1
115               ENDIF
116            END DO
117         END DO
118         !
119         IF(lwp) THEN
120            WRITE(numout,*) '        Namelist nam_diaharm'
121            WRITE(numout,*) '        nb_ana    = ', nb_ana
122            CALL flush(numout)
123         ENDIF
124         !
125         IF (nb_ana > jpmax_harmo) THEN
126            WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo'
127            WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo
128            CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 )
129         ENDIF
130
131         ALLOCATE(name    (nb_ana))
132         DO jh=1,nb_ana
133            DO ji=1,jpmax_harmo
134               IF (TRIM(tname(jh)) ==  Wave(ji)%cname_tide) THEN
135                  name(jh) = ji
136                  EXIT
137               END IF
138            END DO
139         END DO
140
141         ! Initialize frequency array:
142         ! ---------------------------
143         ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
144
145         CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
146
147         IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
148
149         DO jh = 1, nb_ana
150            IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
151         END DO
152
153         ! Initialize temporary arrays:
154         ! ----------------------------
155         ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
156         ana_temp(:,:,:,:) = 0._wp
157
158      ENDIF
159
160   END SUBROUTINE dia_harm_init
161
162
163   SUBROUTINE dia_harm ( kt )
164      !!----------------------------------------------------------------------
165      !!                 ***  ROUTINE dia_harm  ***
166      !!         
167      !! ** Purpose :   Tidal harmonic analysis main routine
168      !!
169      !! ** Action  :   Sums ssh/u/v over time analysis [nit000_han,nitend_han]
170      !!
171      !!--------------------------------------------------------------------
172      INTEGER, INTENT( IN ) :: kt
173      !
174      INTEGER  :: ji, jj, jh, jc, nhc
175      REAL(wp) :: ztime, ztemp
176      !!--------------------------------------------------------------------
177      IF( ln_timing )   CALL timing_start('dia_harm')
178      !
179      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
180         !
181         ztime = (kt-nit000+1) * rdt 
182         !
183         nhc = 0
184         DO jh = 1, nb_ana
185            DO jc = 1, 2
186               nhc = nhc+1
187               ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
188                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
189                  !
190               DO jj = 1,jpj
191                  DO ji = 1,jpi
192                     ! Elevation
193                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)       
194                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
195                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
196                  END DO
197               END DO
198               !
199            END DO
200         END DO
201         !       
202      END IF
203      !
204      IF( kt == nitend_han )   CALL dia_harm_end
205      !
206      IF( ln_timing )   CALL timing_stop('dia_harm')
207      !
208   END SUBROUTINE dia_harm
209
210
211   SUBROUTINE dia_harm_end
212      !!----------------------------------------------------------------------
213      !!                 ***  ROUTINE diaharm_end  ***
214      !!         
215      !! ** Purpose :  Compute the Real and Imaginary part of tidal constituents
216      !!
217      !! ** Action  :  Decompose the signal on the harmonic constituents
218      !!
219      !!--------------------------------------------------------------------
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,jpmax_harmo,2) ::   ana_amp   ! workspace
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      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 )/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               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      ! Elevation:
258      DO jj = 1, jpj
259         DO ji = 1, jpi
260            ! Fill input array
261            kun = 0
262            DO jh = 1, nb_ana
263               DO jc = 1, 2
264                  kun = kun + 1
265                  ztmp4(kun)=ana_temp(ji,jj,kun,1)
266               END DO
267            END DO
268
269            CALL SUR_DETERMINE(jj)
270
271            ! Fill output array
272            DO jh = 1, nb_ana
273               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
274               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
275            END DO
276         END DO
277      END DO
278
279      ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   & 
280         &      out_u  (jpi,jpj,2*nb_ana),   &
281         &      out_v  (jpi,jpj,2*nb_ana)  )
282
283      DO jj = 1, jpj
284         DO ji = 1, jpi
285            DO jh = 1, nb_ana 
286               X1 = ana_amp(ji,jj,jh,1)
287               X2 =-ana_amp(ji,jj,jh,2)
288               out_eta(ji,jj,jh       ) = X1 * tmask_i(ji,jj)
289               out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj)
290            END DO
291         END DO
292      END DO
293
294      ! ubar:
295      DO jj = 1, jpj
296         DO ji = 1, jpi
297            ! Fill input array
298            kun=0
299            DO jh = 1,nb_ana
300               DO jc = 1,2
301                  kun = kun + 1
302                  ztmp4(kun)=ana_temp(ji,jj,kun,2)
303               END DO
304            END DO
305
306            CALL SUR_DETERMINE(jj+1)
307
308            ! Fill output array
309            DO jh = 1, nb_ana
310               ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
311               ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
312            END DO
313
314         END DO
315      END DO
316
317      DO jj = 1, jpj
318         DO ji = 1, jpi
319            DO jh = 1, nb_ana 
320               X1= ana_amp(ji,jj,jh,1)
321               X2=-ana_amp(ji,jj,jh,2)
322               out_u(ji,jj,       jh) = X1 * ssumask(ji,jj)
323               out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj)
324            ENDDO
325         ENDDO
326      ENDDO
327
328      ! vbar:
329      DO jj = 1, jpj
330         DO ji = 1, jpi
331            ! Fill input array
332            kun=0
333            DO jh = 1,nb_ana
334               DO jc = 1,2
335                  kun = kun + 1
336                  ztmp4(kun)=ana_temp(ji,jj,kun,3)
337               END DO
338            END DO
339
340            CALL SUR_DETERMINE(jj+1)
341
342            ! Fill output array
343            DO jh = 1, nb_ana
344               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
345               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
346            END DO
347
348         END DO
349      END DO
350
351      DO jj = 1, jpj
352         DO ji = 1, jpi
353            DO jh = 1, nb_ana 
354               X1=ana_amp(ji,jj,jh,1)
355               X2=-ana_amp(ji,jj,jh,2)
356               out_v(ji,jj,       jh)=X1 * ssvmask(ji,jj)
357               out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj)
358            END DO
359         END DO
360      END DO
361      !
362      CALL dia_wri_harm ! Write results in files
363      !
364   END SUBROUTINE dia_harm_end
365
366
367   SUBROUTINE dia_wri_harm
368      !!--------------------------------------------------------------------
369      !!                 ***  ROUTINE dia_wri_harm  ***
370      !!         
371      !! ** Purpose : Write tidal harmonic analysis results in a netcdf file
372      !!--------------------------------------------------------------------
373      CHARACTER(LEN=lc) :: cltext
374      CHARACTER(LEN=lc) ::   &
375         cdfile_name_T   ,   & ! name of the file created (T-points)
376         cdfile_name_U   ,   & ! name of the file created (U-points)
377         cdfile_name_V         ! name of the file created (V-points)
378      INTEGER  ::   jh
379      !!----------------------------------------------------------------------
380
381      IF(lwp) WRITE(numout,*) '  '
382      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
383      IF(lwp) WRITE(numout,*) '  '
384
385      ! A) Elevation
386      !/////////////
387      !
388      DO jh = 1, nb_ana
389      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
390      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
391      END DO
392
393      ! B) ubar
394      !/////////
395      !
396      DO jh = 1, nb_ana
397      CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
398      CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) )
399      END DO
400
401      ! C) vbar
402      !/////////
403      !
404      DO jh = 1, nb_ana
405         CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) )
406         CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
407      END DO
408      !
409   END SUBROUTINE dia_wri_harm
410
411
412   SUBROUTINE SUR_DETERMINE(init)
413      !!---------------------------------------------------------------------------------
414      !!                      *** ROUTINE SUR_DETERMINE ***
415      !!   
416      !!   
417      !!       
418      !!---------------------------------------------------------------------------------
419      INTEGER, INTENT(in) ::   init 
420      !
421      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd
422      REAL(wp)                        :: zval1, zval2, zx1
423      REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2
424      INTEGER , DIMENSION(jpincomax) :: ipos2, ipivot
425      !---------------------------------------------------------------------------------
426      !           
427      IF( init == 1 ) THEN
428         IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
429         IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
430         !
431         ztmp3(:,:) = 0._wp
432         !
433         DO jh1_sd = 1, nsparse
434            DO jh2_sd = 1, nsparse
435               nisparse(jh2_sd) = nisparse(jh2_sd)
436               njsparse(jh2_sd) = njsparse(jh2_sd)
437               IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN
438                  ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  &
439                     &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd)
440               ENDIF
441            END DO
442         END DO
443         !
444         DO jj_sd = 1 ,ninco
445            ipos1(jj_sd) = jj_sd
446            ipos2(jj_sd) = jj_sd
447         END DO
448         !
449         DO ji_sd = 1 , ninco
450            !
451            !find greatest non-zero pivot:
452            zval1 = ABS(ztmp3(ji_sd,ji_sd))
453            !
454            ipivot(ji_sd) = ji_sd
455            DO jj_sd = ji_sd, ninco
456               zval2 = ABS(ztmp3(ji_sd,jj_sd))
457               IF( zval2 >= zval1 )THEN
458                  ipivot(ji_sd) = jj_sd
459                  zval1         = zval2
460               ENDIF
461            END DO
462            !
463            DO ji1_sd = 1, ninco
464               zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd)
465               zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd))
466               ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd)
467               ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
468            END DO
469            !
470            ipos2(ji_sd)         = ipos1(ipivot(ji_sd))
471            ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
472            ipos1(ji_sd)         = ipos2(ji_sd)
473            ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
474            zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd)
475            DO jj_sd = 1, ninco
476               ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
477            END DO
478            !
479            DO ji2_sd = ji_sd+1, ninco
480               zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
481               DO jj_sd=1,ninco
482                  ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
483               END DO
484            END DO
485            !
486         END DO
487         !
488      ENDIF ! End init==1
489
490      DO ji_sd = 1, ninco
491         ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
492         DO ji2_sd = ji_sd+1, ninco
493            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
494         END DO
495      END DO
496
497      !system solving:
498      ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
499      ji_sd = ninco
500      DO ji_sd = ninco-1, 1, -1
501         zx1 = 0._wp
502         DO jj_sd = ji_sd+1, ninco
503            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
504         END DO
505         ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
506      END DO
507
508      DO jj_sd =1, ninco
509         ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
510      END DO
511      !
512   END SUBROUTINE SUR_DETERMINE
513
514   !!======================================================================
515END MODULE diaharm
Note: See TracBrowser for help on using the repository browser.