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/r14075_India_uncoupled/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm.F90 @ 15422

Last change on this file since 15422 was 15422, checked in by jcastill, 11 months ago

Changes tested so that they can merged with the CO9 Met Office branch - jpmax_harmo should be 34 with FES14 tides, but the last components are not used anyway

File size: 28.8 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   !!   NB: 2017-12 : add 3D harmonic analysis of velocities 
9   !!                 integration of Maria Luneva's development 
10   !!   'key_3Ddiaharm
11   !!----------------------------------------------------------------------
12   USE oce             ! ocean dynamics and tracers variables
13   USE dom_oce         ! ocean space and time domain
14   USE phycst
15   USE daymod
16   USE tide_mod
17   USE sbctide         ! Tidal forcing or not
18   !
19# if defined key_3Ddiaharm 
20   USE zdf_oce 
21#endif 
22   !
23   USE in_out_manager  ! I/O units
24   USE iom             ! I/0 library
25   USE ioipsl          ! NetCDF IPSL library
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE timing          ! preformance summary
28   USE lib_mpp           ! MPP library
29
30   IMPLICIT NONE
31   PRIVATE
32   
33   INTEGER, PARAMETER :: jpincomax    = 2.*jpmax_harmo
34   INTEGER, PARAMETER :: jpdimsparse  = jpincomax*366*24*2 ! 30min for a 1yr-long run
35
36   !                         !!** namelist variables **
37   LOGICAL, PUBLIC ::   ln_diaharm    ! Choose tidal harmonic output or not
38   INTEGER         ::   nit000_han    ! First time step used for harmonic analysis
39   INTEGER         ::   nitend_han    ! Last time step used for harmonic analysis
40   INTEGER         ::   nstep_han     ! Time step frequency for harmonic analysis
41   INTEGER         ::   nb_ana        ! Number of harmonics to analyse
42
43   INTEGER , ALLOCATABLE, DIMENSION(:)           ::   name 
44   REAL(wp), ALLOCATABLE, DIMENSION(:)           ::   ana_freq, ut, vt, ft 
45# if defined key_3Ddiaharm 
46   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:)   ::   ana_temp 
47   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     ::   out_eta, out_u, out_v, out_w, out_dzi 
48# else 
49   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     ::   ana_temp 
50   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       ::   out_eta, out_u, out_v 
51# endif
52
53   INTEGER  ::   ninco, nsparse
54   REAL(wp) ::   z1_tmp3
55   INTEGER ,       DIMENSION(jpdimsparse)         ::   njsparse, nisparse
56   INTEGER , SAVE, DIMENSION(jpincomax)           ::   ipos1
57   REAL(wp),       DIMENSION(jpdimsparse)         ::   valuesparse
58   REAL(wp),       DIMENSION(jpincomax)           ::   ztmp4 , ztmp7, z1_pivot
59   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) ::   ztmp3 , zpilier
60
61   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
62
63   PUBLIC   dia_harm        ! routine called by step.F90
64   PUBLIC   dia_harm_init   ! routine called by nemogcm.F90
65
66   !!----------------------------------------------------------------------
67   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
68   !! $Id$
69   !! Software governed by the CeCILL license (see ./LICENSE)
70   !!----------------------------------------------------------------------
71CONTAINS
72
73   SUBROUTINE dia_harm_init 
74      !!----------------------------------------------------------------------
75      !!                 ***  ROUTINE dia_harm_init  ***
76      !!         
77      !! ** Purpose :   Initialization of tidal harmonic analysis
78      !!
79      !! ** Method  :   Initialize frequency array and  nodal factor for nit000_han
80      !!
81      !!--------------------------------------------------------------------
82      INTEGER ::   jh, nhan, ji
83      INTEGER ::   ios                 ! Local integer output status for namelist read
84
85      NAMELIST/nam_diaharm/ ln_diaharm, nit000_han, nitend_han, nstep_han, tname
86      !!----------------------------------------------------------------------
87
88      IF(lwp) THEN
89         WRITE(numout,*)
90         WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
91# if defined key_3Ddiaharm 
92         WRITE(numout,*) '  - 3D harmonic analysis of currents activated (key_3Ddiaharm)' 
93#endif
94         WRITE(numout,*) '~~~~~~~~~~~~~ '
95      ENDIF
96      !
97      REWIND( numnam_ref )              ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
98      READ  ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
99901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_diaharm in reference namelist' )
100      REWIND( numnam_cfg )              ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis
101      READ  ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 )
102902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist' )
103      IF(lwm) WRITE ( numond, nam_diaharm )
104      !
105      IF(lwp) THEN
106         WRITE(numout,*) 'Tidal diagnostics = ', ln_diaharm
107         WRITE(numout,*) '   First time step used for analysis:         nit000_han= ', nit000_han
108         WRITE(numout,*) '   Last  time step used for analysis:         nitend_han= ', nitend_han
109         WRITE(numout,*) '   Time step frequency for harmonic analysis: nstep_han = ', nstep_han
110      ENDIF
111
112      IF( ln_diaharm .AND. .NOT.ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
113
114      IF( ln_diaharm ) THEN
115
116         CALL tide_init_Wave
117         !
118         ! Basic checks on harmonic analysis time window:
119         ! ----------------------------------------------
120         IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   &
121            &                                       ' restart capability not implemented' )
122         IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   &
123            &                                       'restart capability not implemented' )
124
125         IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   &
126            &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
127         !
128         nb_ana = 0
129         DO jh=1,jpmax_harmo
130            DO ji=1,jpmax_harmo
131               IF(TRIM(tname(jh)) == Wave(ji)%cname_tide) THEN
132                  nb_ana=nb_ana+1
133               ENDIF
134            END DO
135         END DO
136         !
137         IF(lwp) THEN
138            WRITE(numout,*) '        Namelist nam_diaharm'
139            WRITE(numout,*) '        nb_ana    = ', nb_ana
140            CALL flush(numout)
141         ENDIF
142         !
143         IF (nb_ana > jpmax_harmo) THEN
144            WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo'
145            WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo
146            CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 )
147         ENDIF
148
149         ALLOCATE(name(nb_ana))
150         DO jh=1,nb_ana
151            DO ji=1,jpmax_harmo
152               IF (TRIM(tname(jh)) ==  Wave(ji)%cname_tide) THEN
153                  name(jh) = ji
154                  EXIT
155               END IF
156            END DO
157         END DO
158
159         ! Initialize frequency array:
160         ! ---------------------------
161         ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
162
163         CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
164
165         IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
166
167         DO jh = 1, nb_ana
168            IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
169         END DO
170
171         ! Initialize temporary arrays:
172         ! ----------------------------
173# if defined key_3Ddiaharm 
174         ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 5, jpk ) ) 
175         ana_temp(:,:,:,:,:) = 0._wp 
176# else 
177         ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 3      ) ) 
178         ana_temp(:,:,:,:  ) = 0._wp 
179#endif
180
181      ENDIF
182
183   END SUBROUTINE dia_harm_init
184
185
186   SUBROUTINE dia_harm( kt )
187      !!----------------------------------------------------------------------
188      !!                 ***  ROUTINE dia_harm  ***
189      !!         
190      !! ** Purpose :   Tidal harmonic analysis main routine
191      !!
192      !! ** Action  :   Sums ssh/u/v over time analysis [nit000_han,nitend_han]
193      !!
194      !!--------------------------------------------------------------------
195      INTEGER, INTENT( in ) ::   kt
196      !
197      INTEGER  ::   ji, jj, jh, jc, nhc
198# if defined key_3Ddiaharm 
199      INTEGER  :: jk 
200# endif
201      REAL(wp) ::   ztime, ztemp
202      !!--------------------------------------------------------------------
203      IF( ln_timing )   CALL timing_start('dia_harm')
204      !
205      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
206         !
207         ztime = (kt-nit000+1) * rdt 
208         !
209         nhc = 0
210         DO jh = 1, nb_ana
211            DO jc = 1, 2
212               nhc = nhc+1
213               ztemp =  (   MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
214                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
215                  !
216               ! ssh, ub, vb are stored at the last level of 5d array
217               DO jj = 2, jpjm1
218                  DO ji = 2, jpim1
219                     ! Elevation and currents 
220# if defined key_3Ddiaharm 
221                     ana_temp(ji,jj,nhc,1,jpk) = ana_temp(ji,jj,nhc,1,jpk) + ztemp*sshn(ji,jj)*ssmask (ji,jj)         
222                     ana_temp(ji,jj,nhc,2,jpk) = ana_temp(ji,jj,nhc,2,jpk) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 
223                     ana_temp(ji,jj,nhc,3,jpk) = ana_temp(ji,jj,nhc,3,jpk) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 
224 
225                     ana_temp(ji,jj,nhc,5,jpk) = ana_temp(ji,jj,nhc,5,jpk)                               & 
226                   &                              + ztemp*bfrva(ji,jj)*vn(ji,jj,mbkv(ji,jj))*ssvmask(ji,jj) 
227                     ana_temp(ji,jj,nhc,4,jpk) = ana_temp(ji,jj,nhc,4,jpk)                               &   
228                   &                              + ztemp*bfrua(ji,jj)*un(ji,jj,mbku(ji,jj))*ssumask(ji,jj) 
229# else 
230                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)         
231                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 
232                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 
233# endif
234                  END DO
235               END DO
236               !
237# if defined key_3Ddiaharm 
238! 3d velocity and density: 
239             DO jk=1,jpk-1 
240               DO jj = 1,jpj 
241                  DO ji = 1,jpi 
242                     ! density and velocity 
243                     ana_temp(ji,jj,nhc,1,jk) = ana_temp(ji,jj,nhc,1,jk) + ztemp*rhd(ji,jj,jk) 
244                     ana_temp(ji,jj,nhc,2,jk) = ana_temp(ji,jj,nhc,2,jk) + ztemp*(un(ji,jj,jk)-un_b(ji,jj)) & 
245                &                                          *umask(ji,jj,jk) 
246                     ana_temp(ji,jj,nhc,3,jk) = ana_temp(ji,jj,nhc,3,jk) + ztemp*(vn(ji,jj,jk)-vn_b(ji,jj)) & 
247                &                                          *vmask(ji,jj,jk)   
248                     ana_temp(ji,jj,nhc,4,jk) = ana_temp(ji,jj,nhc,4,jk) + ztemp*wn(ji,jj,jk) 
249   
250                     ana_temp(ji,jj,nhc,5,jk) = ana_temp(ji,jj,nhc,5,jk) - 0.5*grav*ztemp*(rhd(ji,jj,jk)+rhd(ji,jj,jk+1))/max(rn2(ji,jj,jk),1.e-8_wp) 
251                  END DO 
252               END DO 
253             ENDDO 
254# endif
255            END DO
256         END DO
257      END IF
258      !
259      IF( kt == nitend_han )   CALL dia_harm_end
260      !
261      IF( ln_timing )   CALL timing_stop('dia_harm')
262      !
263   END SUBROUTINE dia_harm
264
265
266   SUBROUTINE dia_harm_end
267      !!----------------------------------------------------------------------
268      !!                 ***  ROUTINE diaharm_end  ***
269      !!         
270      !! ** Purpose :  Compute the Real and Imaginary part of tidal constituents
271      !!
272      !! ** Action  :  Decompose the signal on the harmonic constituents
273      !!
274      !!--------------------------------------------------------------------
275      INTEGER  ::   ji, jj, jh, jc, jn, nhan
276# if defined key_3Ddiaharm 
277      INTEGER :: jk 
278# endif
279      INTEGER  ::   ksp, kun, keq
280      REAL(wp) ::   ztime, ztime_ini, ztime_end, z1_han
281      !!--------------------------------------------------------------------
282      !
283      IF(lwp) WRITE(numout,*)
284      IF(lwp) WRITE(numout,*) 'dia_harm_end: kt=nitend_han: Perform harmonic analysis'
285      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
286     
287      ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis
288      ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis
289      nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
290      z1_han = 1._wp / REAL(nhan-1) 
291     
292# if defined key_3Ddiaharm 
293      ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana),   & 
294         &      out_u  (jpi,jpj,jpk,2*nb_ana),   & 
295         &      out_v  (jpi,jpj,jpk,2*nb_ana),   & 
296         &      out_w  (jpi,jpj,jpk,2*nb_ana),   & 
297         &      out_dzi(jpi,jpj,jpk,2*nb_ana) ) 
298# else 
299      ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   & 
300         &      out_u  (jpi,jpj,2*nb_ana),   & 
301         &      out_v  (jpi,jpj,2*nb_ana)  ) 
302# endif 
303 
304      IF(lwp) WRITE(numout,*) 'ANA F OLD', ft   
305      IF(lwp) WRITE(numout,*) 'ANA U OLD', ut 
306      IF(lwp) WRITE(numout,*) 'ANA V OLD', vt
307
308      ninco = 2*nb_ana
309
310      ksp = 0
311      keq = 0       
312      DO jn = 1, nhan
313         ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end ) * z1_han
314         keq = keq + 1
315         kun = 0
316         DO jh = 1, nb_ana
317            DO jc = 1, 2
318               kun = kun + 1
319               ksp = ksp + 1
320               nisparse(ksp) = keq
321               njsparse(ksp) = kun
322               valuesparse(ksp) = (   MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))   &
323                  &             + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
324            END DO
325         END DO
326      END DO
327
328      nsparse = ksp
329
330      IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
331      IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
332
333      CALL SUR_DETERMINE_INIT
334
335      ! Density and Elevation: 
336# if defined key_3Ddiaharm 
337    DO jk=1,jpk 
338# endif
339      DO jj = 2, jpjm1
340         DO ji = 2, jpim1
341
342            ! Fill input array
343# if defined key_3Ddiaharm 
344            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1,jk) 
345# else
346            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1)
347# endif
348            CALL SUR_DETERMINE
349           
350            ! Fill output array
351            DO jh = 1, nb_ana
352# if defined key_3Ddiaharm 
353               out_eta(ji,jj,jk,jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj)
354               out_eta(ji,jj,jk,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj)
355# else 
356               out_eta(ji,jj,   jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj)
357               out_eta(ji,jj,   jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj)
358# endif
359            END DO
360         END DO
361      END DO
362
363      ! ubar:
364      DO jj = 2, jpjm1
365         DO ji = 2, jpim1
366
367            ! Fill input array
368# if defined key_3Ddiaharm 
369            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2,jk)
370# else
371            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2)
372# endif
373            CALL SUR_DETERMINE
374
375            ! Fill output array
376            DO jh = 1, nb_ana
377# if defined key_3Ddiaharm 
378               out_u(ji,jj,jk,       jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj)
379               out_u(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj)
380# else 
381               out_u(ji,jj,          jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj)
382               out_u(ji,jj,   nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj)
383# endif
384
385            END DO
386
387        END DO
388      END DO
389
390      ! vbar:
391      DO jj = 2, jpjm1
392         DO ji = 2, jpim1
393
394            ! Fill input array
395# if defined key_3Ddiaharm 
396            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3,jk)
397# else
398            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3)
399# endif
400            CALL SUR_DETERMINE
401
402            ! Fill output array
403            DO jh = 1, nb_ana
404# if defined key_3Ddiaharm 
405               out_v(ji,jj,jk,       jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj)
406               out_v(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj)
407# else 
408               out_v(ji,jj,          jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj)
409               out_v(ji,jj,   nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj)
410# endif
411            END DO
412
413         END DO
414      END DO
415
416# if defined key_3Ddiaharm 
417      ! w- velocity 
418      DO jj = 1, jpj 
419         DO ji = 1, jpi 
420            ! Fill input array 
421            kun=0 
422            DO jh = 1,nb_ana 
423               DO jc = 1,2 
424                  kun = kun + 1 
425                  ztmp4(kun)=ana_temp(ji,jj,kun,4,jk) 
426               END DO 
427            END DO 
428 
429            CALL SUR_DETERMINE(jj+1) 
430 
431            ! Fill output array 
432            DO jh = 1, nb_ana 
433               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 
434               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 
435            END DO 
436 
437         END DO 
438      END DO 
439 
440      DO jj = 1, jpj 
441         DO ji = 1, jpi 
442            DO jh = 1, nb_ana 
443               X1=ana_amp(ji,jj,jh,1) 
444               X2=-ana_amp(ji,jj,jh,2) 
445               out_w(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj) 
446               out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj) 
447            END DO 
448         END DO 
449      END DO 
450 
451       ! dzi- isopycnal displacements 
452      DO jj = 1, jpj 
453         DO ji = 1, jpi 
454            ! Fill input array 
455            kun=0 
456            DO jh = 1,nb_ana 
457               DO jc = 1,2 
458                  kun = kun + 1 
459                  ztmp4(kun)=ana_temp(ji,jj,kun,5,jk) 
460               END DO 
461            END DO 
462 
463            CALL SUR_DETERMINE(jj+1) 
464 
465            ! Fill output array 
466            DO jh = 1, nb_ana 
467               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 
468               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 
469            END DO 
470 
471         END DO 
472      END DO 
473 
474      DO jj = 1, jpj 
475         DO ji = 1, jpi 
476            DO jh = 1, nb_ana 
477               X1=ana_amp(ji,jj,jh,1) 
478               X2=-ana_amp(ji,jj,jh,2) 
479               out_dzi(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj) 
480               out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj) 
481            END DO 
482         END DO 
483      END DO 
484 
485   ENDDO ! jk 
486# endif
487      !
488      ! clem: we could avoid this call if all the loops were from 1:jpi and 1:jpj
489      !       but I think this is the most efficient
490      CALL lbc_lnk_multi( 'dia_harm_end', out_eta, 'T', 1., out_u, 'U', -1. , out_v, 'V', -1. )
491      !
492      CALL dia_wri_harm ! Write results in files
493      !
494      DEALLOCATE( out_eta, out_u, out_v )
495      !
496   END SUBROUTINE dia_harm_end
497
498
499   SUBROUTINE dia_wri_harm
500      !!--------------------------------------------------------------------
501      !!                 ***  ROUTINE dia_wri_harm  ***
502      !!         
503      !! ** Purpose : Write tidal harmonic analysis results in a netcdf file
504      !!--------------------------------------------------------------------
505      INTEGER  ::   jh
506 
507# if defined key_3Ddiaharm 
508      CHARACTER(LEN=lc) :: cdfile_name_W         ! name of the file created (W-points) 
509      INTEGER  :: jk 
510      REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im   
511      REAL(WP), ALLOCATABLE, DIMENSION (:,:)   :: z2real, z2im       
512# endif 
513!!---------------------------------------------------------------------- 
514 
515#if defined key_dimgout 
516      cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc' 
517      cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc' 
518      cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc' 
519#   if defined key_3Ddiaharm 
520      cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc' 
521#   endif 
522#endif
523
524      IF(lwp) WRITE(numout,*) '  '
525      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
526#if defined key_dimgout 
527      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T) 
528      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_U) 
529      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_V) 
530#   if defined key_3Ddiaharm 
531      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_W) 
532#   endif 
533#endif
534      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
535
536# if defined key_3Ddiaharm 
537      ALLOCATE(z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj)) 
538# endif 
539 
540      ! A) density and elevation
541      !/////////////
542#if defined key_dimgout 
543      cltext='density amplitude and phase; elevation is level=jpk ' 
544      CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2') 
545#else 
546#   if defined key_3Ddiaharm 
547      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 
548#   endif
549      DO jh = 1, nb_ana
550#   if defined key_3Ddiaharm 
551        DO jk=1,jpkm1 
552          z3real(:,:,jk)=out_eta(:,:,jk,jh) 
553          z3im  (:,:,jk)=out_eta(:,:,jk,jh+nb_ana) 
554        ENDDO 
555      z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana) 
556      CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) ) 
557      CALL iom_put( TRIM(tname(jh))//'y_ro', z3im  (:,:,:) ) 
558      CALL iom_put( TRIM(tname(jh))//'x'   , z2real(:,:  ) ) 
559      CALL iom_put( TRIM(tname(jh))//'y'   , z2im  (:,:  ) ) 
560#   else   
561      WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh))
562      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
563      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,jh+nb_ana) )
564#   endif 
565      END DO 
566#endif 
567 
568      ! B) u
569      !/////////
570#if defined key_dimgout 
571      cltext='3d u amplitude and phase; ubar is the last level' 
572      CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2') 
573#else 
574#   if defined key_3Ddiaharm 
575      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 
576#   endif
577      DO jh = 1, nb_ana
578#   if defined key_3Ddiaharm 
579        DO jk=1,jpkm1 
580          z3real(:,:,jk)=out_u(:,:,jk,jh) 
581          z3im  (:,:,jk)=out_u(:,:,jk,jh+nb_ana) 
582        ENDDO 
583        z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana) 
584        CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) ) 
585        CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:)  ) 
586        CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) ) 
587        CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:)  ) 
588        z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana) 
589        CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) ) 
590        CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:)  ) 
591#   else 
592        CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) ) 
593        CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) ) 
594#   endif 
595      END DO 
596#endif 
597 
598      ! C) v
599      !/////////
600#if defined key_dimgout 
601      cltext='3d v amplitude and phase; vbar is the last level' 
602      CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2') 
603#else 
604#   if defined key_3Ddiaharm 
605      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 
606#   endif
607      DO jh = 1, nb_ana
608#   if defined key_3Ddiaharm 
609        DO jk=1,jpkm1 
610          z3real(:,:,jk)=out_v(:,:,jk,jh) 
611          z3im  (:,:,jk)=out_v(:,:,jk,jh+nb_ana) 
612        ENDDO 
613        z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana) 
614        CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) ) 
615        CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:)  ) 
616        CALL iom_put( TRIM(tname(jh))//'x_v2d'  , z2real(:,:) ) 
617        CALL iom_put( TRIM(tname(jh))//'y_v2d'  , z2im (:,:)  ) 
618        z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana) 
619        CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) ) 
620        CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:)  ) 
621#   else 
622         CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) ) 
623         CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) ) 
624#   endif 
625       END DO 
626 
627#endif 
628      ! D) w 
629# if defined key_3Ddiaharm 
630#   if defined key_dimgout 
631      cltext='3d w amplitude and phase; vort_baro is the last level' 
632      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2') 
633#   else 
634      DO jh = 1, nb_ana 
635        DO jk=1,jpkm1 
636         z3real(:,:,jk)=out_w(:,:,jk,jh) 
637         z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana) 
638        ENDDO 
639        CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) ) 
640        CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) ) 
641      END DO 
642#   endif 
643 
644!       E) dzi + tau_bot 
645#   if defined key_dimgout 
646      cltext='dzi=g*ro/N2 amplitude and phase' 
647      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2') 
648#   else 
649      DO jh = 1, nb_ana 
650        DO jk=1,jpkm1 
651         z3real(:,:,jk)=out_dzi(:,:,jk,jh) 
652         z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana) 
653        ENDDO 
654        CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) ) 
655        CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) ) 
656      END DO 
657#   endif 
658# endif   
659 
660     
661# if defined key_3Ddiaharm 
662   DEALLOCATE(z3real, z3im, z2real,z2im) 
663# endif
664
665   END SUBROUTINE dia_wri_harm
666
667
668   SUBROUTINE SUR_DETERMINE_INIT
669      !!---------------------------------------------------------------------------------
670      !!                      *** ROUTINE SUR_DETERMINE_INIT ***
671      !!       
672      !!---------------------------------------------------------------------------------
673      INTEGER                        :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd
674      INTEGER                        :: ipivot
675      REAL(wp)                       :: zval1, zval2, zcol1, zcol2
676      INTEGER , DIMENSION(jpincomax) :: ipos2
677      !!---------------------------------------------------------------------------------
678      !           
679      !
680      ztmp3(:,:) = 0._wp
681      !
682      DO jh1_sd = 1, nsparse
683         DO jh2_sd = 1, nsparse
684            IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN
685               ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  &
686                  &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd)
687            ENDIF
688         END DO
689      END DO
690      !
691      DO jj_sd = 1, ninco
692         ipos1(jj_sd) = jj_sd
693         ipos2(jj_sd) = jj_sd
694      END DO
695      !
696      DO ji_sd = 1, ninco
697         !
698         !find greatest non-zero pivot:
699         zval1 = ABS(ztmp3(ji_sd,ji_sd))
700         !
701         ipivot = ji_sd
702         DO jj_sd = ji_sd, ninco
703            zval2 = ABS(ztmp3(ji_sd,jj_sd))
704            IF( zval2 >= zval1 )THEN
705               ipivot = jj_sd
706               zval1  = zval2
707            ENDIF
708         END DO
709         !
710         DO ji1_sd = 1, ninco
711            zcol1                = ztmp3(ji1_sd,ji_sd)
712            zcol2                = ztmp3(ji1_sd,ipivot)
713            ztmp3(ji1_sd,ji_sd)  = zcol2
714            ztmp3(ji1_sd,ipivot) = zcol1
715         END DO
716         !
717         ipos2(ji_sd)    = ipos1(ipivot)
718         ipos2(ipivot)   = ipos1(ji_sd)
719         ipos1(ji_sd)    = ipos2(ji_sd)
720         ipos1(ipivot)   = ipos2(ipivot)
721         z1_pivot(ji_sd) = 1._wp / ztmp3(ji_sd,ji_sd)
722         DO jj_sd = 1, ninco
723            ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) * z1_pivot(ji_sd)
724         END DO
725         !
726         DO ji2_sd = ji_sd+1, ninco
727            zpilier(ji2_sd,ji_sd) = ztmp3(ji2_sd,ji_sd)
728            DO jj_sd=1,ninco
729               ztmp3(ji2_sd,jj_sd) = ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
730            END DO
731         END DO
732         !
733      END DO
734      !
735      z1_tmp3 = 1._wp / ztmp3(ninco,ninco)
736      !
737   END SUBROUTINE SUR_DETERMINE_INIT
738
739   
740   SUBROUTINE SUR_DETERMINE
741      !!---------------------------------------------------------------------------------
742      !!                      *** ROUTINE SUR_DETERMINE ***
743      !!   
744      !!---------------------------------------------------------------------------------
745      INTEGER                        :: ji_sd, jj_sd, ji1_sd, ji2_sd
746      REAL(wp)                       :: zx1
747      REAL(wp), DIMENSION(jpincomax) :: ztmpx
748      !!---------------------------------------------------------------------------------
749      !           
750      DO ji_sd = 1, ninco
751         ztmp4(ji_sd) = ztmp4(ji_sd) * z1_pivot(ji_sd)
752         DO ji2_sd = ji_sd+1, ninco
753            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
754         END DO
755      END DO
756
757      !system solving:
758      ztmpx(ninco) = ztmp4(ninco) * z1_tmp3
759      DO ji_sd = ninco-1, 1, -1
760         zx1 = 0._wp
761         DO jj_sd = ji_sd+1, ninco
762            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
763         END DO
764         ztmpx(ji_sd) = ztmp4(ji_sd) - zx1
765      END DO
766
767      DO jj_sd = 1, ninco
768         ztmp7(ipos1(jj_sd)) = ztmpx(jj_sd)
769      END DO
770      !
771   END SUBROUTINE SUR_DETERMINE
772
773   
774   !!======================================================================
775END MODULE diaharm
Note: See TracBrowser for help on using the repository browser.