source: NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DIA/diaharm_fast.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: 34.9 KB
Line 
1MODULE diaharm_fast 
2   !!======================================================================
3   !!                       ***  MODULE  example  ***
4   !! Ocean physics:  On line harmonic analyser
5   !!                 
6   !!=====================================================================
7
8#if defined key_diaharm_fast 
9
10   !!----------------------------------------------------------------------
11   !!   'key_harm_ana'  :                Calculate harmonic analysis
12   !!----------------------------------------------------------------------
13   !!   harm_ana        :
14   !!   harm_ana_init   :
15   !!   NB: 2017-12 : add 3D harmonic analysis of velocities
16   !!                 integration of Maria Luneva's development
17   !!   'key_3Ddiaharm'
18   !!----------------------------------------------------------------------
19
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE iom
23   USE in_out_manager  ! I/O units
24   USE phycst          ! physical constants
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE bdy_oce         ! ocean open boundary conditions
27   USE bdytides        ! tidal bdy forcing
28   USE zdfdrg  , ONLY : rCdU_bot   ! bottom friction
29   USE daymod          ! calendar
30   USE tideini
31   USE restart
32   USE ioipsl, ONLY : ju2ymds    ! for calendar
33   !
34   !
35   USE timing          ! preformance summary
36   USE zdf_oce
37
38   IMPLICIT NONE
39   PRIVATE
40
41   !! *  Routine accessibility
42   PUBLIC dia_harm_fast                                      ! routine called in step.F90 module
43   LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast  = .TRUE.   ! to be run or not
44   LOGICAL, PUBLIC :: lk_diaharm_2D   ! = .TRUE.   ! to run 2d
45   LOGICAL, PUBLIC :: lk_diaharm_3D   ! = .TRUE.   ! to run 3d
46
47   !! * Module variables
48   INTEGER, PARAMETER ::  nharm_max  = jpmax_harmo  ! max number of harmonics to be analysed
49   INTEGER, PARAMETER ::  nhm_max    = 2*nharm_max+1 
50   INTEGER, PARAMETER ::  nvab       = 2 ! number of 3D variables
51   INTEGER            ::  nharm
52   INTEGER            ::  nhm 
53   INTEGER ::                 & !!! ** toto namelist (namtoto) **
54      nflag  =  1                ! default value of nflag
55   REAL(wp), DIMENSION(nharm_max) ::                & 
56      om_tide                     ! tidal frequencies ( rads/sec)
57   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:)   ::                & 
58      bzz,c,x    ! work arrays
59   REAL(wp) :: cca,ssa,zm,bt,dd_cumul
60!
61   REAL(wp), PUBLIC ::   fjulday_startharm       !: Julian Day since start of harmonic analysis
62   REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf   ! nodel/phase corrections used by diaharmana
63   REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:)   :: cc,a
64!
65   INTEGER ::  nvar_2d, nvar_3d    !: number of 2d and 3d variables to analyse
66   INTEGER, ALLOCATABLE,DIMENSION(:) :: m_posi_2d, m_posi_3d
67
68!  Name of variables used in the restart
69   CHARACTER( LEN = 10 ), DIMENSION(5), PARAMETER :: m_varName2d = (/'ssh','u2d','v2d','ubfr','vbfr'/)
70   CHARACTER( LEN = 10 ), DIMENSION(4), PARAMETER :: m_varName3d = (/'rho','u3d','v3d','w3d'/)
71!
72   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:  ) :: g_cosamp2D, g_sinamp2D, g_cumul_var2D
73   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:) :: g_cosamp3D, g_sinamp3D, g_cumul_var3D 
74!
75   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:)       :: g_out2D,h_out2D  ! arrays for output
76   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:)     :: g_out3D,h_out3D  ! arrays for 3D output
77!
78!  NAMELIST
79   LOGICAL, PUBLIC :: ln_diaharm_store           !: =T  Stores data for harmonic Analysis
80   LOGICAL, PUBLIC :: ln_diaharm_compute         !: =T  Compute harmonic Analysis
81   LOGICAL, PUBLIC :: ln_diaharm_read_restart   !: =T  Read restart from a previous run
82   LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d
83   INTEGER ::   nb_ana        ! Number of harmonics to analyse
84   CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
85   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   ntide_all ! INDEX within the full set of constituents (tide.h90)
86   INTEGER , ALLOCATABLE, DIMENSION(:)       ::   ntide_sub ! INDEX within the subset of constituents pass in input
87
88   !! * Substitutions
89#  include "vectopt_loop_substitute.h90"
90
91   !!----------------------------------------------------------------------
92   !!  OPA 9.0 , LOCEAN-IPSL (2005)
93   !! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
94   !! or  TOP 1.0 , LOCEAN-IPSL (2005)
95   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
96   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
97   !!----------------------------------------------------------------------
98
99CONTAINS
100
101   SUBROUTINE dia_harm_fast( kt )
102      !!----------------------------------------------------------------------
103      !!                    ***  ROUTINE harm_ana  ***
104      !!
105      !! ** Purpose :   Harmonic analyser
106      !!
107      !! ** Method  :   
108      !!
109      !! ** Action  : - first action (share memory array/varible modified
110      !!                in this routine
111      !!              - second action .....
112      !!              - .....
113      !!
114      !! References :
115      !!   Give references if exist otherwise suppress these lines
116      !!
117      !! History :
118      !!   9.0  !  03-08  (Autor Names)  Original code
119      !!        !  02-08  (Author names)  brief description of modifications
120      !!----------------------------------------------------------------------
121      !! * Modules used
122     
123      !! * arguments
124      INTEGER, INTENT( in  ) ::   & 
125         kt                          ! describe it!!!
126
127      !! * local declarations
128      INTEGER  :: ji, jk, jj          ! dummy loop arguments
129      INTEGER  :: jh, i1, i2, jgrid
130      INTEGER  :: j2d, j3d
131      REAL(WP) :: sec2start
132      !!--------------------------------------------------------------------
133
134      IF( ln_timing )   CALL timing_start( 'dia_harm_fast' )
135      IF( kt == nit000   )   CALL harm_ana_init    ! Initialization (first time-step only)
136
137     IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
138
139      ! this bit done every time step
140      nhm=2*nb_ana+1
141      c(1) = 1.0
142
143      sec2start = nint( (fjulday-fjulday_startharm)*86400._wp ) 
144
145      DO jh=1,nb_ana
146         c(2*jh  ) = anaf(jh)*cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
147         c(2*jh+1) = anaf(jh)*sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
148      ENDDO 
149
150      ! CUMULATE
151       DO jj = 2, jpjm1
152          DO ji = fs_2, fs_jpim1   ! vector opt.
153            DO jh=1,nhm   ! loop harmonic
154
155               DO j2d=1,nvar_2d
156                  IF ( m_posi_2d(j2d) .eq. 1 ) dd_cumul = c(jh) * sshn(ji,jj) * ssmask (ji,jj)             ! analysis elevation
157                  IF ( m_posi_2d(j2d) .eq. 2 ) dd_cumul = c(jh) * un_b(ji,jj) * ssumask(ji,jj)             ! analysis depth average velocities
158                  IF ( m_posi_2d(j2d) .eq. 3 ) dd_cumul = c(jh) * vn_b(ji,jj) * ssvmask(ji,jj)
159                  IF ( m_posi_2d(j2d) .eq. 4 ) dd_cumul = c(jh) * 0.5*(rCdU_bot(ji+1,jj) + rCdU_bot(ji,jj)) * &
160                                                          un(ji,jj,mbku(ji,jj)) * ssumask(ji,jj) ! analysis bottom friction
161                  IF ( m_posi_2d(j2d) .eq. 5 ) dd_cumul = c(jh) * 0.5*(rCdU_bot(ji,jj+1) + rCdU_bot(ji,jj)) * &
162                                                          vn(ji,jj,mbkv(ji,jj)) * ssvmask(ji,jj)
163                  g_cumul_var2D(jh,ji,jj,j2d) = g_cumul_var2D(jh,ji,jj,j2d) + dd_cumul
164               ENDDO
165
166               DO j3d=1,nvar_3d
167                  DO jk=1,jpkm1
168                     IF ( m_posi_3d(j3d) .eq. 1 ) dd_cumul = c(jh) *  rhd(ji,jj,jk)               * tmask(ji,jj,jk)   
169                     IF ( m_posi_3d(j3d) .eq. 2 ) dd_cumul = c(jh) * ( un(ji,jj,jk)-un_b(ji,jj) ) * umask(ji,jj,jk) 
170                     IF ( m_posi_3d(j3d) .eq. 3 ) dd_cumul = c(jh) * ( vn(ji,jj,jk)-vn_b(ji,jj) ) * vmask(ji,jj,jk)
171                     IF ( m_posi_3d(j3d) .eq. 4 ) dd_cumul = c(jh) *   wn(ji,jj,jk)               * wmask(ji,jj,jk)
172                     g_cumul_var3D(jh,ji,jj,jk,j3d) = g_cumul_var3D(jh,ji,jj,jk,j3d) + dd_cumul
173                  ENDDO
174               ENDDO
175
176            ENDDO     ! end loop harmonic
177         ENDDO        ! end loop lat
178      ENDDO           ! end loop lon
179
180      ! Compute nodal factor cumulative cross-product
181      DO i1=1,nhm
182         DO i2=1,nhm
183            cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
184         ENDDO
185      ENDDO
186
187      ! Output RESTART
188      IF( kt == nitrst ) THEN
189         CALL harm_rst_write(kt) ! Dump out data for a restarted run
190      ENDIF
191
192      ! At End of run
193      IF ( kt ==  nitend ) THEN
194
195         IF(lwp) WRITE(numout,*)
196         IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
197         IF(lwp) WRITE(numout,*) '~~~~~~~~~'
198
199         IF( ln_diaharm_compute ) THEN
200
201             ! INITIALISE TABLE TO 0
202             IF ( nvar_2d .gt. 0 ) THEN
203                g_cosamp2D = 0.0_wp
204                g_sinamp2D = 0.0_wp
205             ENDIF
206             IF ( nvar_3d .gt. 0 ) THEN
207                g_cosamp3D = 0.0_wp
208                g_sinamp3D = 0.0_wp
209             ENDIF
210
211             ! FIRST OUTPUT 2D VARIABLES
212             DO jgrid=1,nvar_2d    ! loop number of 2d variables (ssh, U2d, V2d, UVfric) to analyse harmonically
213                DO ji=1,jpi        ! loop lon
214                   DO jj=1,jpj     ! loop lat
215                      bt = 1.0_wp; bzz(:) = 0.0_wp
216                      DO jh=1,nhm  ! loop harmonic
217                         bzz(jh) = g_cumul_var2D(jh,ji,jj,jgrid)
218                         bt = bt*bzz(jh)
219                      ENDDO
220                      ! Copy back original cumulated nodal factor
221                      a(:,:) = cc(:,:)
222!                     now do gaussian elimination of the system
223!                     a * x = b
224!                     the matrix x is (a0,a1,b1,a2,b2 ...)
225!                     the matrix a and rhs b solved here for x
226                      x=0.0_wp
227                      IF(bt.ne.0.) THEN
228                        CALL gelim( a, bzz, x, nhm )
229!                       Backup output in variables
230                        DO jh=1,nb_ana
231                           g_cosamp2D(jh,ji,jj,jgrid) = x(jh*2  )
232                           g_sinamp2D(jh,ji,jj,jgrid) = x(jh*2+1)
233                        ENDDO
234                        g_cosamp2D( 0,ji,jj,jgrid) = x(1)
235                        g_sinamp2D( 0,ji,jj,jgrid) = 0.0_wp
236                      ENDIF     ! bt.ne.0.
237                   ENDDO        ! jj
238                ENDDO           ! ji
239             ENDDO              ! jgrid
240
241             ! SECOND OUTPUT 3D VARIABLES
242             DO jgrid=1,nvar_3d     ! loop number of 3d variables rho, U, V, W
243                DO jk=1,jpkm1       ! loop over vertical level
244                   DO ji=1,jpi      ! loop over lon
245                      DO jj=1,jpj   ! loop over lat
246                         bt = 1.0_wp; bzz(:) = 0.0_wp
247                         DO jh=1,nhm
248                            bzz(jh) = g_cumul_var3D(jh,ji,jj,jk,jgrid)
249                            bt = bt*bzz(jh)
250                         ENDDO
251                         ! Copy back original cumulated nodal factor
252                         a(:,:) = cc(:,:)                     
253!                        now do gaussian elimination of the system
254!                        a * x = b
255!                        the matrix x is (a0,a1,b1,a2,b2 ...)
256!                        the matrix a and rhs b solved here for x
257                         x=0.0_wp
258                         IF(bt.ne.0.) THEN
259                           CALL gelim( a, bzz, x, nhm )
260!                          Backup output in variables
261                           DO jh=1,nb_ana
262                              g_cosamp3D(jh,ji,jj,jk,jgrid) = x(jh*2  )
263                              g_sinamp3D(jh,ji,jj,jk,jgrid) = x(jh*2+1)
264                           ENDDO
265                           g_cosamp3D   ( 0,ji,jj,jk,jgrid) = x(1)
266                           g_sinamp3D   ( 0,ji,jj,jk,jgrid) = 0.0_wp
267                        ENDIF     ! bt.ne.0.
268                      ENDDO       ! jj
269                   ENDDO          ! ji
270                ENDDO             ! jk
271             ENDDO                ! jgrid
272
273             CALL harm_ana_out     ! output analysis (last time step)
274
275         ELSE    ! ln_harmana_compute = False
276             IF(lwp) WRITE(numout,*) " Skipping Computing harmonics at last step"
277
278         ENDIF   ! ln_harmana_compute
279      ENDIF      ! kt ==  nitend
280
281     ENDIF
282
283      IF( ln_timing )   CALL timing_stop( 'dia_harm_fast' )
284
285   END SUBROUTINE dia_harm_fast 
286
287   SUBROUTINE harm_ana_init
288      !!----------------------------------------------------------------------
289      !!                  ***  ROUTINE harm_ana_init  ***
290      !!                   
291      !! ** Purpose :   initialization of ....
292      !!
293      !! ** Method  :   blah blah blah ...
294      !!
295      !! ** input   :   Namlist namexa
296      !!
297      !! ** Action  :   ... 
298      !!
299      !! history :
300      !!   9.0  !  03-08  (Autor Names)  Original code
301      !!----------------------------------------------------------------------
302      !! * local declarations
303      INTEGER ::   ji, jk, jh  ! dummy loop indices
304      INTEGER ::   ios                  ! Local integer output status for namelist read
305      INTEGER ::   k2d, k3d             ! dummy number of analysis
306      NAMELIST/nam_diaharm_fast/ ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, tname
307      !!----------------------------------------------------------------------
308
309      lk_diaharm_2D    = .TRUE.   ! to run 2d
310      lk_diaharm_3D    = .TRUE.   ! to run 3d
311
312      IF(lwp) WRITE(numout,*)
313      IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
314      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
315
316      ! GET NAMELIST DETAILS
317      REWIND( numnam_ref )              ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis
318      READ  ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901)
319901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist' )
320
321      REWIND( numnam_cfg )              ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis
322      READ  ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 )
323902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist' )
324      IF(lwm) WRITE ( numond, nam_diaharm_fast )
325
326      ! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90
327      nb_ana = 0
328      DO jk=1,jpmax_harmo
329         DO ji=1,nb_harmo
330            IF(TRIM(tname(jk)) == Wave( ntide(ji) )%cname_tide ) THEN
331               nb_ana=nb_ana+1
332            ENDIF
333         END DO
334      END DO
335      !
336      IF(lwp) THEN
337         WRITE(numout,*) '        Namelist nam_diaharm_fast'
338         WRITE(numout,*) '        nb_ana    = ', nb_ana
339         CALL flush(numout)
340      ENDIF
341      !
342      IF (nb_ana > nharm_max) THEN
343        IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop'
344        IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max
345        nstop = nstop + 1
346      ENDIF
347
348      ALLOCATE(ntide_all(nb_ana))
349      ALLOCATE(ntide_sub(nb_ana))
350
351      DO jk=1,nb_ana
352       DO ji=1,nb_harmo
353          IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN
354             ntide_sub(jk) = ji
355             ntide_all(jk) = ntide(ji)
356             EXIT
357          END IF
358       END DO
359      END DO
360
361      ! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS
362      nvar_2d = 0; nvar_3d = 0
363      IF ( ln_ana_ssh   ) nvar_2d = nvar_2d + 1       ! analysis elevation
364      IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2       ! analysis depth-averaged velocity
365      IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2       ! analysis bottom friction
366           
367      IF ( ln_ana_rho   ) nvar_3d = nvar_3d + 1       ! analysis density
368      IF ( ln_ana_uv3d  ) nvar_3d = nvar_3d + 2       ! analysis 3D horizontal velocities
369      IF ( ln_ana_w3d   ) nvar_3d = nvar_3d + 1       ! analysis 3D vertical velocity
370
371      ! CHECK IF SOMETHING TO RUN
372      IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE.   ! no 2d to run
373      IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE.   ! no 3d to run
374
375      IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
376
377         ! DO ALLOCATIONS
378         IF ( lk_diaharm_2D ) THEN
379            ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj,    nvar_2d) )
380            ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj,    nvar_2d) )
381            ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj,    nvar_2d) )
382            ALLOCATE( g_out2D (jpi,jpj) )
383            ALLOCATE( h_out2D (jpi,jpj) )
384            ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0
385         ENDIF
386 
387         IF ( lk_diaharm_3D ) THEN
388            ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
389            ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
390            ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
391            ALLOCATE( g_out3D (jpi,jpj,jpk) )
392            ALLOCATE( h_out3D (jpi,jpj,jpk) )
393            ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0
394         ENDIF
395
396         ALLOCATE( cc(nb_ana*2+1,nb_ana*2+1) )
397         ALLOCATE( a (nb_ana*2+1,nb_ana*2+1) )
398         ALLOCATE( bzz(nb_ana*2+1) )
399         ALLOCATE( x  (nb_ana*2+1) )
400         ALLOCATE( c  (nb_ana*2+1) )
401         ALLOCATE( anau(nb_ana) )
402         ALLOCATE( anav(nb_ana) )
403         ALLOCATE( anaf(nb_ana) )
404         ! END ALLOCATE
405
406         ! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL
407         ! MAKES THINGS EASIER AND FASTER LATER
408         ! !!! UGLY !!!
409         jh = 1; k2d = 0; 
410         IF ( ln_ana_ssh   ) THEN
411            k2d = k2d + 1; m_posi_2d(k2d) = jh
412            IF(lwp) WRITE(numout,*) "   - ssh harmonic analysis activated (ln_ana_ssh)"
413         ENDIF
414         jh = jh + 1
415         IF ( ln_ana_uvbar ) THEN
416            k2d = k2d + 1; m_posi_2d(k2d) = jh
417            jh  = jh  + 1 
418            k2d = k2d + 1; m_posi_2d(k2d) = jh
419            IF(lwp) WRITE(numout,*) "   - barotropic currents harmonic analysis activated (ln_ana_uvbar)"
420         ELSE
421            jh  = jh  + 1
422         ENDIF
423         jh = jh + 1
424         IF ( ln_ana_bfric ) THEN
425            k2d = k2d + 1; m_posi_2d(k2d) = jh
426            jh  = jh  + 1; 
427            k2d = k2d + 1; m_posi_2d(k2d) = jh
428            IF(lwp) WRITE(numout,*) "   - bottom friction harmonic analysis activated (ln_ana_vbfr)"
429         ELSE
430            jh  = jh  + 1
431         ENDIF
432
433         ! and for 3D
434         jh = 1; k3d = 0; 
435         IF ( ln_ana_rho  ) THEN
436            k3d = k3d + 1; m_posi_3d(k3d) = jh
437            IF(lwp) WRITE(numout,*) "   - 3D density harmonic analysis activated (ln_ana_rho)"
438         ENDIF
439         jh = jh + 1
440         IF ( ln_ana_uv3d )  THEN
441            k3d = k3d + 1; m_posi_3d(k3d) = jh
442            jh  = jh  + 1 
443            k3d = k3d + 1; m_posi_3d(k3d) = jh
444            IF(lwp) WRITE(numout,*) "   - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)"
445         ELSE
446            jh  = jh  + 1
447         ENDIF
448         jh = jh + 1
449         IF ( ln_ana_w3d ) THEN
450            k3d = k3d + 1; m_posi_3d(k3d) = jh
451            IF(lwp) WRITE(numout,*) "   - 3D vertical currents harmonic analysis activated (ln_ana_w3d)"
452         ENDIF
453
454         ! SELECT AND STORE FREQUENCIES
455         IF(lwp)    WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
456         DO jh=1,nb_ana
457            om_tide(jh) = omega_tide( ntide_sub(jh) ) 
458            IF(lwp) WRITE(numout,*) '        - ',tname(jh),' ',om_tide(jh)
459         ENDDO
460
461         ! READ RESTART IF
462         IF ( ln_diaharm_read_restart ) THEN
463            IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run"
464            ! Need to read in  bssh bz, cc anau anav and anaf
465            call harm_rst_read  ! This reads in from the previous day
466                                ! Currrently the data in in assci format
467         ELSE
468
469            IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh "
470 
471            IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,:  ) = 0.0_wp
472            IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp
473            cc           = 0.0_wp
474            a    (:,:)   = 0.0_wp ! NB
475            bzz  (:)     = 0.0_wp
476            x    (:)     = 0.0_wp
477            c    (:)     = 0.0_wp
478            anau (:)     = 0.0_wp
479            anav (:)     = 0.0_wp
480            anaf (:)     = 0.0_wp
481
482            DO jh = 1, nb_ana
483               anau(jh) = utide ( ntide_sub(jh) )
484               anav(jh) = v0tide( ntide_sub(jh) )
485               anaf(jh) = ftide ( ntide_sub(jh) )
486            END DO
487
488            fjulday_startharm=fjulday !Set this at very start and store
489
490            IF (lwp) THEN
491               WRITE(numout,*) '--------------------------'
492               WRITE(numout,*) '   - Output anaf for check'
493               WRITE(numout,*) 'ANA F', anaf
494               WRITE(numout,*) 'ANA U', anau
495               WRITE(numout,*) 'ANA V', anav
496               WRITE(numout,*) fjulday_startharm
497               WRITE(numout,*) '--------------------------'
498            ENDIF
499
500         ENDIF
501
502      ELSE
503
504         IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis"
505
506      ENDIF
507
508   END SUBROUTINE harm_ana_init
509!
510   SUBROUTINE gelim (a,b,x,n)
511      !!----------------------------------------------------------------------
512      !!                    ***  ROUTINE harm_ana  ***
513      !!
514      !! ** Purpose :   Guassian elimination
515      !!
516      !!
517      !! ** Action  : - first action (share memory array/varible modified
518      !!                in this routine
519      !!              - second action .....
520      !!              - .....
521      !!
522      !! References :
523      !!   Give references if exist otherwise suppress these lines
524      !!
525      !! History :
526        implicit none
527!
528        integer  :: n
529        REAL(WP) :: b(nb_ana*2+1), a(nb_ana*2+1,nb_ana*2+1)
530        REAL(WP) :: x(nb_ana*2+1)
531        INTEGER  :: row,col,prow,pivrow,rrow
532        REAL(WP) :: atemp
533        REAL(WP) :: pivot
534        REAL(WP) :: m
535
536        do row=1,n-1
537           pivrow=row
538           pivot=a(row,n-row+1)
539           do prow=row+1,n
540              if (abs(a(prow,n-row+1)).gt.abs(pivot)  ) then
541                 pivot=a(prow,n-row+1)
542                 pivrow=prow
543              endif
544           enddo
545!  swap row and prow
546           if ( pivrow .ne. row ) then
547              atemp=b(pivrow)
548              b(pivrow)=b(row)
549              b(row)=atemp
550              do col=1,n
551                 atemp=a(pivrow,col)
552                 a(pivrow,col)=a(row,col)
553                 a(row,col)=atemp
554              enddo
555           endif
556
557           do rrow=row+1,n
558              if (a(row,row).ne.0) then
559   
560                 m=-a(rrow,n-row+1)/a(row,n-row+1)
561                 do col=1,n
562                    a(rrow,col)=m*a(row,col)+a(rrow,col)
563                 enddo
564                 b(rrow)=m*b(row)+b(rrow)
565              endif
566           enddo
567        enddo
568!  back substitution now
569
570        x(1)=b(n)/a(n,1)
571        do row=n-1,1,-1
572           x(n-row+1)=b(row)
573           do col=1,(n-row)
574              x(n-row+1)=(x(n-row+1)-a(row,col)*x(col)) 
575           enddo
576
577           x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
578        enddo
579
580        return
581   END SUBROUTINE gelim
582
583   SUBROUTINE harm_ana_out
584      !!----------------------------------------------------------------------
585      !!                  ***  ROUTINE harm_ana_init  ***
586      !!                   
587      !! ** Purpose :   initialization of ....
588      !!
589      !! ** Method  :   blah blah blah ...
590      !!
591      !! ** input   :   Namlist namexa
592      !!
593      !! ** Action  :   ... 
594      !!
595      !! history :
596      !!   9.0  !  03-08  (Autor Names)  Original code
597      !!----------------------------------------------------------------------
598        USE dianam          ! build name of file (routine)
599 
600      !! * local declarations
601      INTEGER :: ji, jj, jk, jgrid, jh    ! dummy loop indices
602!      INTEGER :: nh_T
603!      INTEGER :: nid_harm
604!      CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names
605!      CHARACTER (len=40) :: clhstnamu, clhstnamv    ! temporary names
606      CHARACTER (len=40) :: suffix
607!      REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi  ! temporary scalars
608
609      do jgrid=1,nvar_2d
610          do jh=1,nb_ana
611             h_out2D = 0.0
612             g_out2D = 0.0
613             do jj=1,nlcj
614                do ji=1,nlci
615                   cca=g_cosamp2D(jh,ji,jj,jgrid)
616                   ssa=g_sinamp2D(jh,ji,jj,jgrid)
617                   h_out2D(ji,jj)=sqrt(cca**2+ssa**2)
618                   IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
619                      g_out2D(ji,jj)= 0.0_wp
620                   ELSE
621                      g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca)       
622                   ENDIF
623                   IF (h_out2D(ji,jj).ne.0) THEN
624                       h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh)
625                   ENDIF
626                   IF (g_out2D(ji,jj).ne.0) THEN  !Correct and take modulus
627                       g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
628                       if (g_out2D(ji,jj).gt.360.0) then
629                           g_out2D(ji,jj)=g_out2D(ji,jj)-360.0
630                       else if (g_out2D(ji,jj).lt.0.0) then
631                           g_out2D(ji,jj)=g_out2D(ji,jj)+360.0
632                       endif
633                   ENDIF
634                enddo
635             enddo
636             !
637             ! NETCDF OUTPUT
638             suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
639             CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix), h_out2D(:,:) )
640             CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix), g_out2D(:,:) )
641
642          enddo
643      enddo
644!
645! DO THE SAME FOR 3D VARIABLES
646!
647      do jgrid=1,nvar_3d
648          do jh=1,nb_ana
649             h_out3D = 0.0
650             g_out3D = 0.0
651             DO jk=1,jpkm1
652                do jj=1,nlcj
653                   do ji=1,nlci
654                      cca=g_cosamp3D(jh,ji,jj,jk,jgrid)
655                      ssa=g_sinamp3D(jh,ji,jj,jk,jgrid)
656                      h_out3D(ji,jj,jk)=sqrt(cca**2+ssa**2)
657                      IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
658                         g_out3D(ji,jj,jk) = 0.0_wp
659                      ELSE
660                         g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca)
661                      ENDIF
662                      IF (h_out3D(ji,jj,jk).ne.0) THEN
663                          h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh)
664                      ENDIF
665                      IF (g_out3D(ji,jj,jk).ne.0) THEN  !Correct and take modulus
666                          g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
667                          if      (g_out3D(ji,jj,jk).gt.360.0) then
668                                   g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0
669                          else if (g_out3D(ji,jj,jk).lt.0.0) then
670                                   g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0
671                          endif
672                      ENDIF
673                   enddo    ! ji
674                enddo       ! jj
675             ENDDO          ! jk
676             !
677             ! NETCDF OUTPUT
678             suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
679             IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
680             CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix), h_out3D(:,:,:) )
681             CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix), g_out3D(:,:,:) )
682          enddo             ! jh
683      enddo                 ! jgrid
684!
685   END SUBROUTINE harm_ana_out
686!
687   SUBROUTINE harm_rst_write(kt)
688      !!----------------------------------------------------------------------
689      !!                  ***  ROUTINE harm_ana_init  ***
690      !!                   
691      !! ** Purpose :  To write out cummulated Tidal Harmomnic data to file for
692      !!               restarting
693      !!
694      !! ** Method  :   restart files will be dated by default
695      !!
696      !! ** input   :   
697      !!
698      !! ** Action  :   ... 
699      !!
700      !! history :
701      !!   0.0  !  01-16  (Enda O'Dea)  Original code
702      !! ASSUMES  dated file for rose  , can change later to be more generic
703      !!----------------------------------------------------------------------
704      INTEGER, INTENT(in) ::   kt     ! ocean time-step
705      !!
706      INTEGER             ::   jh, j2d, j3d
707      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
708      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
709      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
710      CHARACTER(LEN=250)  ::   clfinal   ! full name
711
712      !restart file
713      DO j2d=1,nvar_2d
714         CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
715         DO jh=1,nb_ana
716            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2  , :, :, j2d ) )
717            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
718         ENDDO
719      ENDDO
720
721      DO j3d=1,nvar_3d
722         CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
723         DO jh=1,nb_ana
724            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2  , :, :, :, j3d ) )
725            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
726         ENDDO
727      ENDDO
728
729      IF(lwp) THEN
730        IF( kt > 999999999 ) THEN ; WRITE(clkt, *       ) kt
731        ELSE                      ; WRITE(clkt, '(i8.8)') kt
732        ENDIF
733        clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
734        clpath = TRIM(cn_ocerst_outdir)
735        IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
736        IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
737
738        WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
739        OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
740        WRITE(66) cc
741        WRITE(66) anau
742        WRITE(66) anav
743        WRITE(66) anaf
744        WRITE(66) fjulday_startharm
745        CLOSE(66)
746        WRITE(numout,*) '----------------------------'
747        WRITE(numout,*) '   harm_rst_write: DONE '
748        WRITE(numout,*) cc
749        WRITE(numout,*) anaf
750        WRITE(numout,*) fjulday_startharm
751        WRITE(numout,*) '----------------------------'
752      ENDIF
753 
754   END SUBROUTINE harm_rst_write
755
756   SUBROUTINE harm_rst_read
757      !!----------------------------------------------------------------------
758      !!                  ***  ROUTINE harm_ana_init  ***
759      !!                   
760      !! ** Purpose :  To read in  cummulated Tidal Harmomnic data to file for
761      !!               restarting
762      !!
763      !! ** Method  :   
764      !!
765      !! ** input   :   
766      !!
767      !! ** Action  :   ... 
768      !!
769      !! history :
770      !!   0.0  !  01-16  (Enda O'Dea)  Original code
771      !! ASSUMES  dated file for rose  , can change later to be more generic
772      !!----------------------------------------------------------------------
773      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
774      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
775      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
776      CHARACTER(LEN=250)  ::   clfinal   ! full name
777      INTEGER             ::   jh, j2d, j3d
778
779      IF( nit000 > 999999999 ) THEN ; WRITE(clkt, *       ) nit000-1
780      ELSE                      ; WRITE(clkt, '(i8.8)') nit000-1
781      ENDIF
782      clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
783      clpath = TRIM(cn_ocerst_outdir)
784      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
785
786      IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
787
788      DO j2d=1,nvar_2d
789         CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
790         IF(lwp) WRITE(numout,*) "2D", j2d, m_posi_2d(j2d), m_varName2d( m_posi_2d(j2d) )
791         DO jh=1,nb_ana
792            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2  , :, :, j2d ) )
793            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
794         ENDDO
795      ENDDO
796
797      DO j3d=1,nvar_3d
798         CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
799         IF(lwp) WRITE(numout,*) "3D", j3d,  m_posi_3d(j3d), m_varName3d( m_posi_3d(j3d) )
800
801         DO jh=1,nb_ana
802            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2  , :, :, :, j3d ) )
803            CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
804         ENDDO
805      ENDDO
806
807      WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
808      OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
809      READ(66) cc
810      READ(66) anau
811      READ(66) anav
812      READ(66) anaf
813      READ(66) fjulday_startharm
814      CLOSE(66)
815
816      IF(lwp) THEN
817        WRITE(numout,*) '----------------------------'
818        WRITE(numout,*) '   Checking anaf is correct'
819        WRITE(numout,*) cc
820        WRITE(numout,*) anaf
821        WRITE(numout,*) fjulday_startharm
822        WRITE(numout,*) '----------------------------'
823      ENDIF
824 
825   END SUBROUTINE harm_rst_read
826
827   !!======================================================================
828#else
829!!---------------------------------------------------------------------------------
830!!   Dummy module                                   NO harmonic Analysis
831!!---------------------------------------------------------------------------------
832        LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast  = .FALSE.   ! to be run or not
833
834        CONTAINS
835           SUBROUTINE harm_rst_write(kt)     ! Dummy routine
836           END SUBROUTINE harm_rst_write
837           SUBROUTINE harm_rst_read    ! Dummy routine
838           END SUBROUTINE harm_rst_read
839           SUBROUTINE harm_ana_out      ! Dummy routine
840           END SUBROUTINE harm_ana_out
841           SUBROUTINE harm_ana_init
842           END SUBROUTINE harm_ana_init
843           SUBROUTINE harm_ana( kt )
844!--- NB : end call not properly written
845           END SUBROUTINE harm_ana
846!           END SUBROUTINE harm_ana_init
847!--- END NB
848           SUBROUTINE gelim (a,b,x,n)
849!--- NB : end call not properly written
850           END SUBROUTINE gelim
851!           END SUBROUTINE gelim (a,b,x,n)
852!--- END NB           
853#endif
854
855END MODULE diaharm_fast 
Note: See TracBrowser for help on using the repository browser.