1 | MODULE 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 | |
---|
44 | CONTAINS |
---|
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 | |
---|
209 | END MODULE diacfl |
---|