source: NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DIA/diaharm.F90 @ 12453

Last change on this file since 12453 was 12453, checked in by jcastill, 12 months ago

First implementation of the branch - compiling after merge

File size: 29.5 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*300*24
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   INTEGER ,       DIMENSION(jpdimsparse)         ::   njsparse, nisparse
55   INTEGER , SAVE, DIMENSION(jpincomax)           ::   ipos1
56   REAL(wp),       DIMENSION(jpdimsparse)         ::   valuesparse
57   REAL(wp),       DIMENSION(jpincomax)           ::   ztmp4 , ztmp7
58   REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) ::   ztmp3 , zpilier
59   REAL(wp), SAVE, DIMENSION(jpincomax)           ::   zpivot
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 = 1,jpj
218                  DO ji = 1,jpi
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         !       
258      END IF
259      !
260      IF( kt == nitend_han )   CALL dia_harm_end
261      !
262      IF( ln_timing )   CALL timing_stop('dia_harm')
263      !
264   END SUBROUTINE dia_harm
265
266
267   SUBROUTINE dia_harm_end
268      !!----------------------------------------------------------------------
269      !!                 ***  ROUTINE diaharm_end  ***
270      !!         
271      !! ** Purpose :  Compute the Real and Imaginary part of tidal constituents
272      !!
273      !! ** Action  :  Decompose the signal on the harmonic constituents
274      !!
275      !!--------------------------------------------------------------------
276      INTEGER :: ji, jj, jh, jc, jn, nhan, jl
277# if defined key_3Ddiaharm 
278      INTEGER :: jk 
279# endif
280      INTEGER :: ksp, kun, keq
281      REAL(wp) :: ztime, ztime_ini, ztime_end
282      REAL(wp) :: X1, X2
283      REAL(wp), DIMENSION(jpi,jpj,jpmax_harmo,2) ::   ana_amp   ! workspace
284      !!--------------------------------------------------------------------
285      !
286      IF(lwp) WRITE(numout,*)
287      IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
288      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
289
290      ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis
291      ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis
292      nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
293
294# if defined key_3Ddiaharm 
295      ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana),   & 
296         &      out_u  (jpi,jpj,jpk,2*nb_ana),   & 
297         &      out_v  (jpi,jpj,jpk,2*nb_ana),   & 
298         &      out_w  (jpi,jpj,jpk,2*nb_ana),   & 
299         &      out_dzi(jpi,jpj,jpk,2*nb_ana) ) 
300# else
301      ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   & 
302         &      out_u  (jpi,jpj,2*nb_ana),   & 
303         &      out_v  (jpi,jpj,2*nb_ana)  ) 
304# endif
305 
306      IF(lwp) WRITE(numout,*) 'ANA F OLD', ft 
307      IF(lwp) WRITE(numout,*) 'ANA U OLD', ut 
308      IF(lwp) WRITE(numout,*) 'ANA V OLD', vt
309
310      ninco = 2*nb_ana
311
312      ksp = 0
313      keq = 0       
314      DO jn = 1, nhan
315         ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end )/FLOAT(nhan-1)
316         keq = keq + 1
317         kun = 0
318         DO jh = 1, nb_ana
319            DO jc = 1, 2
320               kun = kun + 1
321               ksp = ksp + 1
322               nisparse(ksp) = keq
323               njsparse(ksp) = kun
324               valuesparse(ksp) = (   MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))   &
325                  &             + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
326            END DO
327         END DO
328      END DO
329
330      nsparse = ksp
331
332      ! Density and Elevation:
333# if defined key_3Ddiaharm 
334    DO jk=1,jpk 
335# endif
336      DO jj = 1, jpj
337         DO ji = 1, jpi
338            ! Fill input array
339            kun = 0
340            DO jh = 1, nb_ana
341               DO jc = 1, 2
342                  kun = kun + 1
343# if defined key_3Ddiaharm 
344                  ztmp4(kun)=ana_temp(ji,jj,kun,1,jk) 
345# else
346                  ztmp4(kun)=ana_temp(ji,jj,kun,1)
347# endif
348               END DO
349            END DO
350
351            CALL SUR_DETERMINE(jj)
352
353            ! Fill output array
354            DO jh = 1, nb_ana
355               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
356               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
357            END DO
358         END DO
359      END DO
360
361      DO jj = 1, jpj
362         DO ji = 1, jpi
363            DO jh = 1, nb_ana 
364               X1 = ana_amp(ji,jj,jh,1)
365               X2 =-ana_amp(ji,jj,jh,2)
366# if defined key_3Ddiaharm 
367               out_eta(ji,jj,jk,jh       ) = X1 * tmask_i(ji,jj) 
368               out_eta(ji,jj,jk,jh+nb_ana) = X2 * tmask_i(ji,jj) 
369# else
370               out_eta(ji,jj   ,jh       ) = X1 * tmask_i(ji,jj) 
371               out_eta(ji,jj   ,jh+nb_ana) = X2 * tmask_i(ji,jj) 
372# endif
373            END DO
374         END DO
375      END DO 
376 
377      ! u-component of velocity
378      DO jj = 1, jpj
379         DO ji = 1, jpi
380            ! Fill input array
381            kun=0
382            DO jh = 1,nb_ana
383               DO jc = 1,2
384                  kun = kun + 1
385# if defined key_3Ddiaharm 
386                  ztmp4(kun)=ana_temp(ji,jj,kun,2,jk) 
387# else
388                  ztmp4(kun)=ana_temp(ji,jj,kun,2)
389# endif
390               END DO
391            END DO
392
393            CALL SUR_DETERMINE(jj+1)
394
395            ! Fill output array
396            DO jh = 1, nb_ana
397               ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
398               ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
399            END DO
400
401         END DO
402      END DO
403
404      DO jj = 1, jpj
405         DO ji = 1, jpi
406            DO jh = 1, nb_ana 
407               X1= ana_amp(ji,jj,jh,1)
408               X2=-ana_amp(ji,jj,jh,2)
409# if defined key_3Ddiaharm 
410               out_u(ji,jj,jk,       jh) = X1 * ssumask(ji,jj) 
411               out_u(ji,jj,jk,nb_ana+jh) = X2 * ssumask(ji,jj) 
412# else
413               out_u(ji,jj,          jh) = X1 * ssumask(ji,jj) 
414               out_u(ji,jj,   nb_ana+jh) = X2 * ssumask(ji,jj) 
415# endif
416            ENDDO
417         ENDDO
418      ENDDO
419
420      ! v- velocity
421      DO jj = 1, jpj
422         DO ji = 1, jpi
423            ! Fill input array
424            kun=0
425            DO jh = 1,nb_ana
426               DO jc = 1,2
427                  kun = kun + 1
428# if defined key_3Ddiaharm 
429                  ztmp4(kun)=ana_temp(ji,jj,kun,3,jk) 
430# else
431                  ztmp4(kun)=ana_temp(ji,jj,kun,3)
432# endif
433               END DO
434            END DO
435
436            CALL SUR_DETERMINE(jj+1)
437
438            ! Fill output array
439            DO jh = 1, nb_ana
440               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
441               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
442            END DO
443
444         END DO
445      END DO
446
447      DO jj = 1, jpj
448         DO ji = 1, jpi
449            DO jh = 1, nb_ana 
450               X1=ana_amp(ji,jj,jh,1)
451               X2=-ana_amp(ji,jj,jh,2)
452# if defined key_3Ddiaharm 
453               out_v(ji,jj,jk,       jh)=X1 * ssvmask(ji,jj) 
454               out_v(ji,jj,jk,nb_ana+jh)=X2 * ssvmask(ji,jj) 
455# else
456               out_v(ji,jj,          jh)=X1 * ssvmask(ji,jj) 
457               out_v(ji,jj,   nb_ana+jh)=X2 * ssvmask(ji,jj) 
458# endif
459            END DO
460         END DO
461      END DO 
462 
463# if defined key_3Ddiaharm 
464      ! w- velocity
465      DO jj = 1, jpj 
466         DO ji = 1, jpi 
467            ! Fill input array
468            kun=0 
469            DO jh = 1,nb_ana 
470               DO jc = 1,2 
471                  kun = kun + 1 
472                  ztmp4(kun)=ana_temp(ji,jj,kun,4,jk) 
473               END DO
474            END DO
475 
476            CALL SUR_DETERMINE(jj+1) 
477 
478            ! Fill output array
479            DO jh = 1, nb_ana 
480               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 
481               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 
482            END DO
483 
484         END DO
485      END DO
486 
487      DO jj = 1, jpj 
488         DO ji = 1, jpi 
489            DO jh = 1, nb_ana 
490               X1=ana_amp(ji,jj,jh,1) 
491               X2=-ana_amp(ji,jj,jh,2) 
492               out_w(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj) 
493               out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj) 
494            END DO
495         END DO
496      END DO 
497 
498       ! dzi- isopycnal displacements
499      DO jj = 1, jpj 
500         DO ji = 1, jpi 
501            ! Fill input array
502            kun=0 
503            DO jh = 1,nb_ana 
504               DO jc = 1,2 
505                  kun = kun + 1 
506                  ztmp4(kun)=ana_temp(ji,jj,kun,5,jk) 
507               END DO
508            END DO
509 
510            CALL SUR_DETERMINE(jj+1) 
511 
512            ! Fill output array
513            DO jh = 1, nb_ana 
514               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 
515               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 
516            END DO
517 
518         END DO
519      END DO
520 
521      DO jj = 1, jpj 
522         DO ji = 1, jpi 
523            DO jh = 1, nb_ana 
524               X1=ana_amp(ji,jj,jh,1) 
525               X2=-ana_amp(ji,jj,jh,2) 
526               out_dzi(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj) 
527               out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj) 
528            END DO
529         END DO
530      END DO
531 
532   ENDDO ! jk
533# endif
534      !
535      CALL dia_wri_harm ! Write results in files
536      !
537   END SUBROUTINE dia_harm_end
538
539
540   SUBROUTINE dia_wri_harm
541      !!--------------------------------------------------------------------
542      !!                 ***  ROUTINE dia_wri_harm  ***
543      !!         
544      !! ** Purpose : Write tidal harmonic analysis results in a netcdf file
545      !!--------------------------------------------------------------------
546      CHARACTER(LEN=lc) :: cltext
547      CHARACTER(LEN=lc) ::   &
548         cdfile_name_T   ,   & ! name of the file created (T-points)
549         cdfile_name_U   ,   & ! name of the file created (U-points)
550         cdfile_name_V         ! name of the file created (V-points)
551      INTEGER  ::   jh
552
553# if defined key_3Ddiaharm 
554      CHARACTER(LEN=lc) :: cdfile_name_W         ! name of the file created (W-points)
555      INTEGER  :: jk 
556      REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im 
557      REAL(WP), ALLOCATABLE, DIMENSION (:,:)   :: z2real, z2im       
558# endif 
559!!----------------------------------------------------------------------
560 
561#if defined key_dimgout 
562      cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc' 
563      cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc' 
564      cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc' 
565#   if defined key_3Ddiaharm 
566      cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc' 
567#   endif 
568#endif
569
570      IF(lwp) WRITE(numout,*) '  '
571      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
572#if defined key_dimgout 
573      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T) 
574      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_U) 
575      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_V) 
576#   if defined key_3Ddiaharm 
577      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_W) 
578#   endif 
579#endif
580      IF(lwp) WRITE(numout,*) '  '
581
582# if defined key_3Ddiaharm 
583      ALLOCATE(z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj)) 
584# endif 
585 
586      ! A) density and elevation
587      !/////////////
588      !
589#if defined key_dimgout 
590      cltext='density amplitude and phase; elevation is level=jpk ' 
591      CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2') 
592#else 
593#   if defined key_3Ddiaharm 
594      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 
595#   endif
596      DO jh = 1, nb_ana
597#   if defined key_3Ddiaharm 
598        DO jk=1,jpkm1 
599          z3real(:,:,jk)=out_eta(:,:,jk,jh) 
600          z3im  (:,:,jk)=out_eta(:,:,jk,jh+nb_ana) 
601        ENDDO 
602      z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana) 
603      CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) ) 
604      CALL iom_put( TRIM(tname(jh))//'y_ro', z3im  (:,:,:) ) 
605      CALL iom_put( TRIM(tname(jh))//'x'   , z2real(:,:  ) ) 
606      CALL iom_put( TRIM(tname(jh))//'y'   , z2im  (:,:  ) ) 
607#   else 
608      WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh))
609      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
610      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
611#   endif
612      END DO 
613#endif 
614 
615      ! B) u
616      !/////////
617      !
618#if defined key_dimgout 
619      cltext='3d u amplitude and phase; ubar is the last level' 
620      CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2') 
621#else 
622#   if defined key_3Ddiaharm 
623      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 
624#   endif
625      DO jh = 1, nb_ana
626#   if defined key_3Ddiaharm 
627        DO jk=1,jpkm1 
628          z3real(:,:,jk)=out_u(:,:,jk,jh) 
629          z3im  (:,:,jk)=out_u(:,:,jk,jh+nb_ana) 
630        ENDDO 
631        z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana) 
632        CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) ) 
633        CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:)  ) 
634        CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) ) 
635        CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:)  ) 
636        z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana) 
637        CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) ) 
638        CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:)  ) 
639#   else
640        CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) ) 
641        CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) ) 
642#   endif
643      END DO 
644#endif 
645 
646      ! C) v
647      !/////////
648      !
649#if defined key_dimgout 
650      cltext='3d v amplitude and phase; vbar is the last level' 
651      CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2') 
652#else 
653#   if defined key_3Ddiaharm 
654      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 
655#   endif
656      DO jh = 1, nb_ana
657#   if defined key_3Ddiaharm 
658        DO jk=1,jpkm1 
659          z3real(:,:,jk)=out_v(:,:,jk,jh) 
660          z3im  (:,:,jk)=out_v(:,:,jk,jh+nb_ana) 
661        ENDDO 
662        z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana) 
663        CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) ) 
664        CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:)  ) 
665        CALL iom_put( TRIM(tname(jh))//'x_v2d'  , z2real(:,:) ) 
666        CALL iom_put( TRIM(tname(jh))//'y_v2d'  , z2im (:,:)  ) 
667        z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana) 
668        CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) ) 
669        CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:)  ) 
670#   else
671         CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) ) 
672         CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) ) 
673#   endif
674       END DO 
675 
676#endif 
677      ! D) w
678# if defined key_3Ddiaharm 
679#   if defined key_dimgout 
680      cltext='3d w amplitude and phase; vort_baro is the last level' 
681      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2') 
682#   else
683      DO jh = 1, nb_ana 
684        DO jk=1,jpkm1 
685         z3real(:,:,jk)=out_w(:,:,jk,jh) 
686         z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana) 
687        ENDDO 
688        CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) ) 
689        CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) ) 
690      END DO 
691#   endif 
692 
693!       E) dzi + tau_bot
694#   if defined key_dimgout 
695      cltext='dzi=g*ro/N2 amplitude and phase' 
696      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2') 
697#   else
698      DO jh = 1, nb_ana 
699        DO jk=1,jpkm1 
700         z3real(:,:,jk)=out_dzi(:,:,jk,jh) 
701         z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana) 
702        ENDDO 
703        CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) ) 
704        CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) ) 
705      END DO 
706#   endif 
707# endif 
708 
709      !
710# if defined key_3Ddiaharm 
711   DEALLOCATE(z3real, z3im, z2real,z2im) 
712# endif
713
714   END SUBROUTINE dia_wri_harm
715
716
717   SUBROUTINE SUR_DETERMINE(init)
718      !!---------------------------------------------------------------------------------
719      !!                      *** ROUTINE SUR_DETERMINE ***
720      !!   
721      !!   
722      !!       
723      !!---------------------------------------------------------------------------------
724      INTEGER, INTENT(in) ::   init 
725      !
726      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd
727      REAL(wp)                        :: zval1, zval2, zx1
728      REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2
729      INTEGER , DIMENSION(jpincomax) :: ipos2, ipivot
730      !---------------------------------------------------------------------------------
731      !           
732      IF( init == 1 ) THEN
733         IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
734         IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
735         !
736         ztmp3(:,:) = 0._wp
737         !
738         DO jh1_sd = 1, nsparse
739            DO jh2_sd = 1, nsparse
740               nisparse(jh2_sd) = nisparse(jh2_sd)
741               njsparse(jh2_sd) = njsparse(jh2_sd)
742               IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN
743                  ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  &
744                     &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd)
745               ENDIF
746            END DO
747         END DO
748         !
749         DO jj_sd = 1 ,ninco
750            ipos1(jj_sd) = jj_sd
751            ipos2(jj_sd) = jj_sd
752         END DO
753         !
754         DO ji_sd = 1 , ninco
755            !
756            !find greatest non-zero pivot:
757            zval1 = ABS(ztmp3(ji_sd,ji_sd))
758            !
759            ipivot(ji_sd) = ji_sd
760            DO jj_sd = ji_sd, ninco
761               zval2 = ABS(ztmp3(ji_sd,jj_sd))
762               IF( zval2 >= zval1 )THEN
763                  ipivot(ji_sd) = jj_sd
764                  zval1         = zval2
765               ENDIF
766            END DO
767            !
768            DO ji1_sd = 1, ninco
769               zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd)
770               zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd))
771               ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd)
772               ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
773            END DO
774            !
775            ipos2(ji_sd)         = ipos1(ipivot(ji_sd))
776            ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
777            ipos1(ji_sd)         = ipos2(ji_sd)
778            ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
779            zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd)
780            DO jj_sd = 1, ninco
781               ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
782            END DO
783            !
784            DO ji2_sd = ji_sd+1, ninco
785               zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
786               DO jj_sd=1,ninco
787                  ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
788               END DO
789            END DO
790            !
791         END DO
792         !
793      ENDIF ! End init==1
794
795      DO ji_sd = 1, ninco
796         ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
797         DO ji2_sd = ji_sd+1, ninco
798            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
799         END DO
800      END DO
801
802      !system solving:
803      ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
804      ji_sd = ninco
805      DO ji_sd = ninco-1, 1, -1
806         zx1 = 0._wp
807         DO jj_sd = ji_sd+1, ninco
808            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
809         END DO
810         ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
811      END DO
812
813      DO jj_sd =1, ninco
814         ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
815      END DO
816      !
817   END SUBROUTINE SUR_DETERMINE
818
819   !!======================================================================
820END MODULE diaharm
Note: See TracBrowser for help on using the repository browser.