source: branches/2017/dev_r8126_ROBUST08_no_ghost/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90 @ 8979

Last change on this file since 8979 was 8979, checked in by acc, 4 years ago

Branch 2017/dev_r8126_ROBUST08_no_ghost. Merge in trunk changes up to rev 8864 in preparation for the merge. Sette tests OK except for AGRIF restartability but uncertain about the trunks status for this?

  • Property svn:keywords set to Id
File size: 9.1 KB
Line 
1MODULE diacfl
2   !!==============================================================================
3   !!                       ***  MODULE  diacfl  ***
4   !! Output CFL diagnostics to ascii file
5   !!==============================================================================
6   !! History :  1.0  !  2010-03  (E. Blockley)  Original code
7   !!                 !  2014-06  (T Graham) Removed CPP key & Updated to vn3.6
8   !!
9   !!----------------------------------------------------------------------
10   !!   dia_cfl        : Compute and output Courant numbers at each timestep
11   !!----------------------------------------------------------------------
12   USE oce             ! ocean dynamics and active tracers
13   USE dom_oce         ! ocean space and time domain
14   USE lib_mpp         ! distribued memory computing
15   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
16   USE in_out_manager  ! I/O manager
17   USE domvvl     
18   USE timing          ! Performance output
19
20   IMPLICIT NONE
21   PRIVATE
22
23   REAL(wp) :: cu_max, cv_max, cw_max                      ! Run max U Courant number
24   INTEGER, DIMENSION(3) :: cu_loc, cv_loc, cw_loc         ! Run max locations
25   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcu_cfl           ! Courant number arrays
26   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcv_cfl           ! Courant number arrays
27   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcw_cfl           ! Courant number arrays
28
29   INTEGER  :: numcfl                                       ! outfile unit
30   CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"      ! ascii filename
31
32   PUBLIC   dia_cfl       ! routine called by step.F90
33   PUBLIC   dia_cfl_init  ! routine called by nemogcm
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
39   !! $Id$
40   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42
43
44CONTAINS
45
46
47   SUBROUTINE dia_cfl ( kt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dia_cfl  ***
50      !!
51      !! ** Purpose :  Compute the Courant numbers Cu=u*dt/dx and Cv=v*dt/dy
52      !!               and output to ascii file 'cfl_diagnostics.ascii'
53      !!----------------------------------------------------------------------
54
55      INTEGER, INTENT(in) ::  kt                            ! ocean time-step index
56
57      REAL(wp) :: zcu_max, zcv_max, zcw_max                 ! max Courant numbers per timestep
58      INTEGER, DIMENSION(3) :: zcu_loc, zcv_loc, zcw_loc    ! max Courant number locations
59
60      REAL(wp) :: dt                                        ! temporary scalars
61      INTEGER, DIMENSION(3) :: zlocu, zlocv, zlocw          ! temporary arrays
62      INTEGER  :: ji, jj, jk                                ! dummy loop indices
63
64     
65      IF( nn_diacfl == 1) THEN
66         IF( nn_timing == 1 )   CALL timing_start('dia_cfl')
67         ! setup timestep multiplier to account for initial Eulerian timestep
68         IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    dt = rdt
69         ELSE                                        ;    dt = rdt * 2.0
70         ENDIF
71
72             ! calculate Courant numbers
73         DO jk = 1, jpk
74            DO jj = 1, jpj
75               DO ji = 1, fs_jpim1   ! vector opt.
76
77                  ! Courant number for x-direction (zonal current)
78                  zcu_cfl(ji,jj,jk) = ABS(un(ji,jj,jk))*dt/e1u(ji,jj)
79
80                  ! Courant number for y-direction (meridional current)
81                  zcv_cfl(ji,jj,jk) = ABS(vn(ji,jj,jk))*dt/e2v(ji,jj)
82
83                  ! Courant number for z-direction (vertical current)
84                  zcw_cfl(ji,jj,jk) = ABS(wn(ji,jj,jk))*dt/e3w_n(ji,jj,jk)
85               END DO
86            END DO         
87         END DO
88
89         ! calculate maximum values and locations
90         IF( lk_mpp ) THEN
91            CALL mpp_maxloc(zcu_cfl,umask,zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3))
92            CALL mpp_maxloc(zcv_cfl,vmask,zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3))
93            CALL mpp_maxloc(zcw_cfl,tmask,zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3))
94         ELSE
95            zlocu = MAXLOC( ABS( zcu_cfl(:,:,:) ) )
96            zcu_loc(1) = zlocu(1) + nimpp - 1
97            zcu_loc(2) = zlocu(2) + njmpp - 1
98            zcu_loc(3) = zlocu(3)
99            zcu_max = zcu_cfl(zcu_loc(1),zcu_loc(2),zcu_loc(3))
100
101            zlocv = MAXLOC( ABS( zcv_cfl(:,:,:) ) )
102            zcv_loc(1) = zlocv(1) + nimpp - 1
103            zcv_loc(2) = zlocv(2) + njmpp - 1
104            zcv_loc(3) = zlocv(3)
105            zcv_max = zcv_cfl(zcv_loc(1),zcv_loc(2),zcv_loc(3))
106
107            zlocw = MAXLOC( ABS( zcw_cfl(:,:,:) ) )
108            zcw_loc(1) = zlocw(1) + nimpp - 1
109            zcw_loc(2) = zlocw(2) + njmpp - 1
110            zcw_loc(3) = zlocw(3)
111            zcw_max = zcw_cfl(zcw_loc(1),zcw_loc(2),zcw_loc(3))
112         ENDIF
113     
114         ! write out to file
115         IF( lwp ) THEN
116            WRITE(numcfl,FMT='(2x,i4,5x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3)
117            WRITE(numcfl,FMT='(11x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cv', zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3)
118            WRITE(numcfl,FMT='(11x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cw', zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3)
119         ENDIF
120
121         ! update maximum Courant numbers from whole run if applicable
122         IF( zcu_max > cu_max ) THEN
123            cu_max = zcu_max
124            cu_loc = zcu_loc
125         ENDIF
126         IF( zcv_max > cv_max ) THEN
127            cv_max = zcv_max
128            cv_loc = zcv_loc
129         ENDIF
130         IF( zcw_max > cw_max ) THEN
131            cw_max = zcw_max
132            cw_loc = zcw_loc
133         ENDIF
134
135         ! at end of run output max Cu and Cv and close ascii file
136         IF( kt == nitend .AND. lwp ) THEN
137            ! to ascii file
138            WRITE(numcfl,*) '******************************************'
139            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', cu_max, cu_loc(1), cu_loc(2), cu_loc(3)
140            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
141            WRITE(numcfl,*) '******************************************'
142            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', cv_max, cv_loc(1), cv_loc(2), cv_loc(3)
143            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
144            WRITE(numcfl,*) '******************************************'
145            WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', cw_max, cw_loc(1), cw_loc(2), cw_loc(3)
146            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
147            CLOSE( numcfl ) 
148
149            ! to ocean output
150            WRITE(numout,*)
151            WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:'
152            WRITE(numout,*) '~~~~~~~~~~~~'
153            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cu', cu_max, 'at (i, j, k) =   &
154                          &   (', cu_loc(1), cu_loc(2), cu_loc(3), ')'
155            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
156            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cv', cv_max, 'at (i, j, k) =   &
157                          &   (', cv_loc(1), cv_loc(2), cv_loc(3), ')'
158            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
159            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cw', cw_max, 'at (i, j, k) =   &
160                          &   (', cw_loc(1), cw_loc(2), cw_loc(3), ')'
161            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
162
163         ENDIF
164
165         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl')
166      ENDIF
167
168   END SUBROUTINE dia_cfl
169
170   SUBROUTINE dia_cfl_init
171      !!----------------------------------------------------------------------
172      !!                  ***  ROUTINE dia_cfl_init  ***
173      !!                   
174      !! ** Purpose :   create output file, initialise arrays
175      !!----------------------------------------------------------------------
176
177
178      IF( nn_diacfl == 1 ) THEN
179         IF( nn_timing == 1 )   CALL timing_start('dia_cfl_init')
180
181         cu_max=0.0
182         cv_max=0.0
183         cw_max=0.0
184
185         ALLOCATE( zcu_cfl(jpi, jpj, jpk), zcv_cfl(jpi, jpj, jpk), zcw_cfl(jpi, jpj, jpk) )
186
187         zcu_cfl(:,:,:)=0.0
188         zcv_cfl(:,:,:)=0.0
189         zcw_cfl(:,:,:)=0.0
190
191         IF( lwp ) THEN
192            WRITE(numout,*)
193            WRITE(numout,*) 'dia_cfl     : Outputting CFL diagnostics to '//TRIM(clname)
194            WRITE(numout,*) '~~~~~~~~~~~~'
195            WRITE(numout,*)
196
197            ! create output ascii file
198            CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
199            WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k'
200            WRITE(numcfl,*) '******************************************'
201         ENDIF
202
203         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl_init')
204
205      ENDIF
206
207   END SUBROUTINE dia_cfl_init
208
209END MODULE diacfl
Note: See TracBrowser for help on using the repository browser.