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.
diaharmana.F90 in branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diaharmana.F90 @ 8059

Last change on this file since 8059 was 8059, checked in by jgraham, 7 years ago

Merged branches required for AMM15 simulations, see ticket #1904.
Merged branches include:
branches/UKMO/CO6_KD490
branches/UKMO/CO6_Restartable_Tidal_Analysis
branches/UKMO/AMM15_v3_6_STABLE

File size: 23.7 KB
Line 
1MODULE harmonic_analysis
2   !!======================================================================
3   !!                       ***  MODULE  example  ***
4   !! Ocean physics:  On line harmonic analyser
5   !!                 
6   !!=====================================================================
7
8#if defined key_harm_ana
9   !!----------------------------------------------------------------------
10   !!   'key_harm_ana'  :                Calculate harmonic analysis
11   !!----------------------------------------------------------------------
12   !!   harm_ana        :
13   !!   harm_ana_init   :
14   !!----------------------------------------------------------------------
15
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE iom
19   USE in_out_manager  ! I/O units
20   USE phycst          ! physical constants
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE bdy_par         ! Unstructured boundary parameters
23   USE bdy_oce         ! ocean open boundary conditions
24   USE bdytides        ! tidal bdy forcing
25   USE daymod          ! calendar
26   USE tideini
27   USE restart
28   USE ioipsl, ONLY : ju2ymds    ! for calendar
29
30   IMPLICIT NONE
31   PRIVATE
32
33   !! *  Routine accessibility
34   PUBLIC harm_ana    ! routine called in step.F90 module
35
36   !! * Module variables
37   INTEGER, PARAMETER ::  nharm_max  =   jptides_max ! max number of harmonics to be analysed
38   INTEGER, PARAMETER ::  nhm_max    =   2*jptides_max+1 
39   INTEGER, PARAMETER ::  nvab       = 2 ! number of 3D variables
40   INTEGER            ::  nharm
41   INTEGER            ::  nhm 
42   INTEGER ::                 & !!! ** toto namelist (namtoto) **
43      nflag  =  1                ! default value of nflag
44   REAL(wp), DIMENSION(nharm_max) ::                & 
45      om_tide                     ! tidal frequencies ( rads/sec)
46   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:)   ::                & 
47      bzz,c,x    ! work arrays
48   REAL(wp) :: cca,ssa,zm,bt
49   REAL(wp) :: ccau,ccav,ssau,ssav
50   REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf   ! nodel/phase corrections used by diaharmana
51!
52   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) ::   &
53      bssh,bubar,bvbar        ! work array for ssh anaylsis
54   REAL(wp), ALLOCATABLE,SAVE, DIMENSION(:,:,:) ::   &
55      cosampz, cosampu,cosampv,      &
56      sinampz, sinampu,sinampv
57   REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a
58!
59   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) ::   &
60      gout,hout        ! arrays for output
61   REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) ::   &
62      huout,hvout,guout,gvout        ! arrays for output
63   REAL(wp), PUBLIC ::   fjulday_startharm       !: Julian Day since start of harmonic analysis
64
65   !! * Substitutions
66
67   !!----------------------------------------------------------------------
68   !!  OPA 9.0 , LOCEAN-IPSL (2005)
69   !! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
70   !! or  TOP 1.0 , LOCEAN-IPSL (2005)
71   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
72   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
73   !!----------------------------------------------------------------------
74
75CONTAINS
76
77   SUBROUTINE harm_ana( kt )
78      !!----------------------------------------------------------------------
79      !!                    ***  ROUTINE harm_ana  ***
80      !!
81      !! ** Purpose :   Harmonic analyser
82      !!
83      !! ** Method  :   
84      !!
85      !! ** Action  : - first action (share memory array/varible modified
86      !!                in this routine
87      !!              - second action .....
88      !!              - .....
89      !!
90      !! References :
91      !!   Give references if exist otherwise suppress these lines
92      !!
93      !! History :
94      !!   9.0  !  03-08  (Autor Names)  Original code
95      !!        !  02-08  (Author names)  brief description of modifications
96      !!----------------------------------------------------------------------
97      !! * Modules used
98     
99      !! * arguments
100      INTEGER, INTENT( in  ) ::   & 
101         kt                          ! describe it!!!
102
103      !! * local declarations
104      INTEGER ::   ji, jj, jk,  &        ! dummy loop arguments
105                   ih,i1,i2,iv,jgrid
106
107
108      !!--------------------------------------------------------------------
109
110
111
112      IF( kt == nit000  )   CALL harm_ana_init    ! Initialization (first time-step only)
113
114! this bit done every time step
115        !nharm=ncpt_anal
116        nharm=nb_harmo
117        nhm=2*nb_harmo+1
118       
119        c(1)=1.0
120
121        do ih=1,nb_harmo
122             c(2*ih)=    cos((fjulday-fjulday_startharm)*86400._wp*om_tide(ih)  )
123             c(2*ih+1)=  sin((fjulday-fjulday_startharm)*86400._wp*om_tide(ih)  )
124        enddo
125
126
127         do ji=1,jpi
128          do jj=1,jpj
129            do ih=1,nhm
130            bssh (ih,ji,jj)=bssh (ih,ji,jj)+c(ih)*sshn(ji,jj)
131            bubar(ih,ji,jj)=bubar(ih,ji,jj)+c(ih)*un_b(ji,jj)
132            bvbar(ih,ji,jj)=bvbar(ih,ji,jj)+c(ih)*vn_b(ji,jj)
133            enddo
134          enddo
135        enddo
136
137!        IF(lwp) WRITE(666,'(4E15.5)') adatrj*86400._wp, sshn(10,10),bssh(2,10,10),bssh(3,10,10)
138!
139         do i1=1,nhm
140           do i2=1,nhm
141             cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
142            enddo
143         enddo
144!
145! At End of run
146
147      IF (kt ==  nitend  ) then
148
149      IF(lwp) WRITE(numout,*)
150      IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
151      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
152      CALL harm_rst_write(kt)     ! Dump out data for a restarted run
153
154
155      if(ln_harm_ana_compute) then
156               WRITE(numout,*) "Computing harmonics at last step"
157!
158        cosampu=0.0_wp
159        sinampu=0.0_wp
160        cosampv=0.0_wp
161        sinampv=0.0_wp
162        cosampz=0.0_wp
163        sinampz=0.0_wp
164!
165     do jgrid=1,3 ! elevation, Ubar,Vbar
166        do ji=1,jpi
167            do jj=1,jpj
168               bt=1.0
169               do ih=1,nhm
170                   if(jgrid .eq. 1) then
171                     bzz(ih)=bssh(ih,ji,jj)
172                   endif
173                   if(jgrid .eq. 2) then
174                     bzz(ih)=bubar(ih,ji,jj)
175                   endif
176                   if(jgrid .eq. 3) then
177                     bzz(ih)=bvbar(ih,ji,jj)
178                   endif
179                   bt=bt*bzz(ih)
180               enddo
181
182               do i1=1,nhm
183                   do i2=1,nhm
184                       a(i1,i2)=cc(i1,i2)
185                   enddo
186               enddo
187!  now do gaussian elimination of the system
188!  a * x = b
189!  the matrix x is (a0,a1,b1,a2,b2 ...)
190!  the matrix a and rhs b solved here for x
191              x=0.0d0
192              if(bt.ne.0.) then
193                  if(lwp .and. ji == 10 .and. ji == 10) then
194              endif
195              call gelim(a,bzz,x,nhm)
196
197              do ih=1,nb_harmo
198                if(jgrid .eq. 1) then
199                   cosampz(ih,ji,jj)=x(ih*2)
200                   sinampz(ih,ji,jj)=x(ih*2+1)
201                endif
202                if(jgrid .eq. 2) then
203                   cosampu(ih,ji,jj)=x(ih*2)
204                   sinampu(ih,ji,jj)=x(ih*2+1)
205                endif
206                if(jgrid .eq. 3) then
207                   cosampv(ih,ji,jj)=x(ih*2)
208                   sinampv(ih,ji,jj)=x(ih*2+1)
209                endif
210              enddo
211                if(jgrid .eq. 1) then
212                  cosampz(0,ji,jj)=x(1)
213                  sinampz(0,ji,jj)=0.0_wp
214                endif
215                if(jgrid .eq. 2) then
216                  cosampu(0,ji,jj)=x(1)
217                  sinampu(0,ji,jj)=0.0_wp
218                endif
219                if(jgrid .eq. 3) then
220                  cosampv(0,ji,jj)=x(1)
221                  sinampv(0,ji,jj)=0.0_wp
222                endif
223              endif
224           enddo
225        enddo
226      enddo ! jgrid
227!
228!
229        CALL harm_ana_out     ! output analysis (last time step)
230     ELSE ! ln_harmana_compute
231         WRITE(numout,*) "Skipping Computing harmonics at last step"
232     ENDIF
233       
234    ENDIF
235   END SUBROUTINE harm_ana
236
237   SUBROUTINE harm_ana_init
238      !!----------------------------------------------------------------------
239      !!                  ***  ROUTINE harm_ana_init  ***
240      !!                   
241      !! ** Purpose :   initialization of ....
242      !!
243      !! ** Method  :   blah blah blah ...
244      !!
245      !! ** input   :   Namlist namexa
246      !!
247      !! ** Action  :   ... 
248      !!
249      !! history :
250      !!   9.0  !  03-08  (Autor Names)  Original code
251      !!----------------------------------------------------------------------
252      !! * local declarations
253      INTEGER ::   ji, jj, jk, jit,ih   ! dummy loop indices
254
255      !!----------------------------------------------------------------------
256      IF(lwp) WRITE(numout,*)
257      IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
258      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
259     ! DO ALLOCATIONS
260
261      ALLOCATE( bubar(nb_harmo*2+1,jpi,jpj) )
262      ALLOCATE( bvbar(nb_harmo*2+1,jpi,jpj) )
263      ALLOCATE( bssh (nb_harmo*2+1,jpi,jpj) )
264
265      ALLOCATE( cosampz(0:nb_harmo*2+1,jpi,jpj))
266      ALLOCATE( cosampu(0:nb_harmo*2+1,jpi,jpj))
267      ALLOCATE( cosampv(0:nb_harmo*2+1,jpi,jpj))
268
269      ALLOCATE( sinampz(0:nb_harmo*2+1,jpi,jpj))
270      ALLOCATE( sinampu(0:nb_harmo*2+1,jpi,jpj))
271      ALLOCATE( sinampv(0:nb_harmo*2+1,jpi,jpj))
272
273      ALLOCATE( cc(nb_harmo*2+1,nb_harmo*2+1) )
274
275      ALLOCATE( a(nb_harmo*2+1,nb_harmo*2+1) )
276
277      ALLOCATE( anau(nb_harmo) )
278      ALLOCATE( anav(nb_harmo) )
279      ALLOCATE( anaf(nb_harmo) )
280
281      ALLOCATE( bzz(nb_harmo*2+1) )
282      ALLOCATE( x(nb_harmo*2+1) )
283      ALLOCATE( c(nb_harmo*2+1) )
284
285      ALLOCATE( gout(jpi,jpj))
286      ALLOCATE( hout(jpi,jpj))
287
288      ALLOCATE( guout(jpi,jpj))
289      ALLOCATE( huout(jpi,jpj))
290
291      ALLOCATE( gvout(jpi,jpj))
292      ALLOCATE( hvout(jpi,jpj))
293
294        ! END ALLOCATE
295
296                 do ih=1,nb_harmo
297                  om_tide(ih)= omega_tide(ih) 
298                 enddo
299       if(ln_harmana_read) then
300               WRITE(numout,*) "Reading previous harmonic data from previous run"
301               ! Need to read in  bssh bz, cc anau anav and anaf
302               call harm_rst_read  ! This reads in from the previous day
303                                   ! Currrently the data in in assci format
304       else
305               WRITE(numout,*) "Starting harmonic analysis from Fresh "
306         bssh(:,:,:)  = 0.0_wp
307         bubar(:,:,:) = 0.0_wp
308         bvbar(:,:,:) = 0.0_wp
309
310         cc=0.0_wp
311
312         anau(:) = 0.0_wp
313         anav(:) = 0.0_wp
314         anaf(:) = 0.0_wp
315         bzz(:) = 0.0_wp
316         x(:) = 0.0_wp
317         c(:) = 0.0_wp
318
319         DO ih = 1, nb_harmo
320            anau(ih) = utide(ih)
321            anav(ih) = v0tide(ih)
322            anaf(ih) = ftide(ih)
323
324         END DO
325         fjulday_startharm=fjulday !Set this at very start and store
326      endif
327
328   END SUBROUTINE harm_ana_init
329!
330   SUBROUTINE gelim (a,b,x,n)
331      !!----------------------------------------------------------------------
332      !!                    ***  ROUTINE harm_ana  ***
333      !!
334      !! ** Purpose :   Guassian elimination
335      !!
336      !! ** Method  :   
337      !!
338      !! ** Action  : - first action (share memory array/varible modified
339      !!                in this routine
340      !!              - second action .....
341      !!              - .....
342      !!
343      !! References :
344      !!   Give references if exist otherwise suppress these lines
345      !!
346      !! History :
347        implicit none
348!
349        integer  :: n
350   REAL(WP) :: b(nb_harmo*2+1),a(nb_harmo*2+1,nb_harmo*2+1)
351   REAL(WP) :: x(nb_harmo*2+1)
352        INTEGER  :: row,col,prow,pivrow,rrow,ntemp
353   REAL(WP) ::  atemp
354   REAL(WP) ::  pivot
355   REAL(WP) ::  m
356   REAL(WP) ::  tempsol(nb_harmo*2+1)
357
358
359        do row=1,n-1
360           pivrow=row
361           pivot=a(row,n-row+1)
362           do prow=row+1,n
363              if (abs(a(prow,n-row+1)).gt.abs(pivot)  ) then
364                 pivot=a(prow,n-row+1)
365                 pivrow=prow
366              endif
367           enddo
368!  swap row and prow
369           if ( pivrow .ne. row ) then
370              atemp=b(pivrow)
371              b(pivrow)=b(row)
372              b(row)=atemp
373              do col=1,n
374                 atemp=a(pivrow,col)
375                 a(pivrow,col)=a(row,col)
376                 a(row,col)=atemp
377              enddo
378           endif
379
380           do rrow=row+1,n
381              if (a(row,row).ne.0) then
382   
383                 m=-a(rrow,n-row+1)/a(row,n-row+1)
384                 do col=1,n
385                    a(rrow,col)=m*a(row,col)+a(rrow,col)
386                 enddo
387                 b(rrow)=m*b(row)+b(rrow)
388              endif
389           enddo
390        enddo
391!  back substitution now
392
393        x(1)=b(n)/a(n,1)
394        do row=n-1,1,-1
395           x(n-row+1)=b(row)
396           do col=1,(n-row)
397              x(n-row+1)=(x(n-row+1)-a(row,col)*x(col)) 
398           enddo
399
400           x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
401        enddo
402
403        return
404        END SUBROUTINE gelim
405
406   SUBROUTINE harm_ana_out
407      !!----------------------------------------------------------------------
408      !!                  ***  ROUTINE harm_ana_init  ***
409      !!                   
410      !! ** Purpose :   initialization of ....
411      !!
412      !! ** Method  :   blah blah blah ...
413      !!
414      !! ** input   :   Namlist namexa
415      !!
416      !! ** Action  :   ... 
417      !!
418      !! history :
419      !!   9.0  !  03-08  (Autor Names)  Original code
420      !!----------------------------------------------------------------------
421        USE dianam          ! build name of file (routine)
422 
423      !! * local declarations
424      INTEGER ::   ji, jj, jk, jgrid,jit, ih, jjk   ! dummy loop indices
425      INTEGER :: nh_T
426      INTEGER :: nid_harm
427      CHARACTER (len=40) ::           &
428         clhstnamt, clop1, clop2          ! temporary names
429      CHARACTER (len=40) ::           &
430         clhstnamu, clhstnamv   ! temporary names
431      REAL(wp) ::   &
432         zsto1, zsto2, zout, zmax,            &  ! temporary scalars
433         zjulian, zdt, zmdi 
434         
435        do jgrid=1,3
436          do ih=1,nb_harmo
437            hout = 0.0
438            gout = 0.0
439            do jj=1,nlcj
440                do ji=1,nlci
441                    if(jgrid .eq. 1) then ! SSH
442                       cca=cosampz(ih,ji,jj)
443                       ssa=sinampz(ih,ji,jj)
444                    endif
445                    if(jgrid .eq. 2) then ! UBAR
446                       cca=cosampu(ih,ji,jj)
447                       ssa=sinampu(ih,ji,jj)
448                    endif
449                    if(jgrid .eq. 3) then ! VBAR
450                       cca=cosampv(ih,ji,jj)
451                       ssa=sinampv(ih,ji,jj)
452                    endif
453
454                    hout(ji,jj)=sqrt(cca**2+ssa**2)
455       
456                    if (cca.ne.0.0) then
457                         gout(ji,jj)=(180.0/rpi)*atan(ssa/cca)
458                    else
459                        if (ssa.ne.0.0) then
460                            if (ssa.gt. 0.0) then
461                                 gout(ji,jj)=90.0
462                            else
463                                 gout(ji,jj)=270.0
464                            endif
465                        else
466                            gout(ji,jj)=0.0
467                        endif
468                    endif
469
470                    if (gout(ji,jj).lt.0) then
471                        gout(ji,jj)=gout(ji,jj)+180.0
472                    endif
473                    if (ssa.lt.0) then
474                        gout(ji,jj)=gout(ji,jj)+180.0
475                    endif
476                   
477                    if (hout(ji,jj).ne.0) then
478                        hout(ji,jj)=hout(ji,jj)/anaf(ih)
479                    endif
480                    if (gout(ji,jj).ne.0) then
481                                                !Convert Rad to degree and take
482                                                !modulus
483                         gout(ji,jj)=gout(ji,jj)+MOD( (anau(ih)+anav(ih))/rad , 360.0)
484                        if (gout(ji,jj).gt.360.0) then
485                            gout(ji,jj)=gout(ji,jj)-360.0
486                        else if (gout(ji,jj).lt.0) then
487                            gout(ji,jj)=gout(ji,jj)+360.0
488                        endif
489                    endif
490                enddo
491            enddo
492!NETCDF OUTPUT
493            if(jgrid==1) then
494             WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'amp'
495             WRITE(numout,*) TRIM(Wave(ntide(ih))%cname_tide)//'phase'
496             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp', hout(:,:) )
497             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase', gout(:,:) )
498            endif
499            if(jgrid==2) then
500             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_u', hout(:,:) )
501             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_u', gout(:,:) )
502            endif
503            if(jgrid==3) then
504             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'amp_v', hout(:,:) )
505             CALL iom_put( TRIM(Wave(ntide(ih))%cname_tide)//'phase_v', gout(:,:) )
506            endif
507          enddo
508        enddo
509 
510!
511   END SUBROUTINE harm_ana_out
512
513   SUBROUTINE harm_rst_write(kt)
514      !!----------------------------------------------------------------------
515      !!                  ***  ROUTINE harm_ana_init  ***
516      !!                   
517      !! ** Purpose :  To write out cummulated Tidal Harmomnic data to file for
518      !!               restarting
519      !!
520      !! ** Method  :   restart files will be dated by default
521      !!
522      !! ** input   :   
523      !!
524      !! ** Action  :   ... 
525      !!
526      !! history :
527      !!   0.0  !  01-16  (Enda O'Dea)  Original code
528      !! ASSUMES  dated file for rose  , can change later to be more generic
529      !!----------------------------------------------------------------------
530      INTEGER, INTENT(in) ::   kt     ! ocean time-step
531      INTEGER             ::   iyear, imonth, iday,ih
532      REAL (wp)           ::   zsec
533      !!
534      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step define as a character
535      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
536      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
537      CHARACTER(LEN=250)  ::   clfinal   ! full name
538
539      CALL ju2ymds( fjulday , iyear, imonth, iday, zsec)
540
541      WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday
542      ! create the file
543      clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana"
544      clpath = TRIM(cn_ocerst_outdir)
545      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
546      WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
547
548
549
550
551!!      clhstnam=TRIM(cexper)//'.restart_harm_ana'
552      ! Open a file to write the data into
553         write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc
554         open(66,file=TRIM(clfinal))
555         write(66,'(1e20.9)') cc
556         !These Three which contain the most data can be moved to a regular
557         !restart file
558
559            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bssh'     , bssh(1,:,:)       )
560            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bubar'    , bubar(1,:,:)       )
561            CALL iom_rstput( kt, nitrst, numrow, 'Mean_bvbar'    , bvbar(1,:,:)       )
562         do ih=1,nb_harmo
563            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos'    , bssh (ih*2  , : , : )     )
564            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin'    , bssh (ih*2+1, : , : )     )
565
566            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos'   , bubar(ih*2  , : , : )    )
567            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin'   , bubar(ih*2+1, : , : )    )
568
569            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos'   , bvbar(ih*2  , : , : )    )
570            CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin'   , bvbar(ih*2+1, : , : )    )
571         enddo
572
573         write(66,'(1e20.9)') anau
574         write(66,'(1e20.9)') anav
575         write(66,'(1e20.9)') anaf
576         write(66,'(1e25.20)') fjulday_startharm
577         close(66)
578 
579   END SUBROUTINE harm_rst_write
580
581   SUBROUTINE harm_rst_read
582      !!----------------------------------------------------------------------
583      !!                  ***  ROUTINE harm_ana_init  ***
584      !!                   
585      !! ** Purpose :  To read in  cummulated Tidal Harmomnic data to file for
586      !!               restarting
587      !!
588      !! ** Method  :   
589      !!
590      !! ** input   :   
591      !!
592      !! ** Action  :   ... 
593      !!
594      !! history :
595      !!   0.0  !  01-16  (Enda O'Dea)  Original code
596      !! ASSUMES  dated file for rose  , can change later to be more generic
597      !!----------------------------------------------------------------------
598      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name
599      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file
600      CHARACTER(LEN=250)  ::   clfinal   ! full name
601      INTEGER             ::   ih
602
603      ! create the file
604      clname = "restart_harm_ana"
605      clpath = TRIM(cn_ocerst_indir)
606      IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
607      WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
608
609
610
611
612!!      clhstnam=TRIM(cexper)//'.restart_harm_ana'
613      ! Open a file to read the data ifrom
614         write (clfinal,'(a,''_'',i4.4)') trim(clpath)//trim(clname),nproc
615         open(66,file=TRIM(clfinal),status='old')
616         read(66,'(1e20.9)') cc
617!2D regular fields can be read from normal restart this saves space and handy to
618!view in netcdf format also.
619         CALL iom_get( numror,jpdom_autoglo, 'Mean_bssh'     , bssh(1,:,:)       )
620         CALL iom_get( numror,jpdom_autoglo, 'Mean_bubar'    , bubar(1,:,:)       )
621         CALL iom_get( numror,jpdom_autoglo, 'Mean_bvbar'    , bvbar(1,:,:)       )
622         
623         do ih=1,nb_harmo
624
625            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_cos'    , bssh (ih*2  , : , : )     )
626            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bssh_sin'    , bssh (ih*2+1, : , : )     )
627
628            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_cos'   , bubar(ih*2  , : , : )    )
629            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bubar_sin'   , bubar(ih*2+1, : , : )    )
630
631            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_cos'   , bvbar(ih*2  , : , : )    )
632            CALL iom_get(   numror,jpdom_autoglo, TRIM(Wave(ntide(ih))%cname_tide)//'bvbar_sin'   , bvbar(ih*2+1, : , : )    )
633         enddo
634         read(66,'(1e20.9)') anau
635         read(66,'(1e20.9)') anav
636         read(66,'(1e20.9)') anaf
637         read(66,'(1e25.20)') fjulday_startharm
638      WRITE(numout,*) 'Checking anaf is correct'
639      WRITE(numout,*) anaf
640      WRITE(numout,*) '---end of anaf checking----'
641
642         close(66)
643
644   END SUBROUTINE harm_rst_read
645
646   !!======================================================================
647#else
648!!---------------------------------------------------------------------------------
649!!   Dummy module                                   NO harmonic Analysis
650!!---------------------------------------------------------------------------------
651        CONTAINS
652           SUBROUTINE harm_rst_write(kt)     ! Dummy routine
653           END SUBROUTINE harm_rst_write
654           SUBROUTINE harm_rst_read    ! Dummy routine
655           END SUBROUTINE harm_rst_read
656           SUBROUTINE harm_ana_out      ! Dummy routine
657           END SUBROUTINE harm_ana_out
658           SUBROUTINE harm_ana_init
659           END SUBROUTINE harm_ana_init
660           SUBROUTINE harm_ana( kt )
661           END SUBROUTINE harm_ana_init
662           SUBROUTINE gelim (a,b,x,n)
663           END SUBROUTINE gelim (a,b,x,n)
664           
665#endif
666
667END MODULE harmonic_analysis
Note: See TracBrowser for help on using the repository browser.