source: branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90 @ 4725

Last change on this file since 4725 was 4725, checked in by timgraham, 6 years ago

Removed unnecessary LBC_LNK calls.

File size: 9.0 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 "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
40   !! $Id$
41   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43
44
45CONTAINS
46
47
48   SUBROUTINE dia_cfl ( kt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE dia_cfl  ***
51      !!
52      !! ** Purpose :  Compute the Courant numbers Cu=u*dt/dx and Cv=v*dt/dy
53      !!               and output to ascii file 'cfl_diagnostics.ascii'
54      !!----------------------------------------------------------------------
55
56      INTEGER, INTENT(in) ::  kt                            ! ocean time-step index
57
58      REAL(wp) :: zcu_max, zcv_max, zcw_max                 ! max Courant numbers per timestep
59      INTEGER, DIMENSION(3) :: zcu_loc, zcv_loc, zcw_loc    ! max Courant number locations
60
61      REAL(wp) :: dt                                        ! temporary scalars
62      INTEGER, DIMENSION(3) :: zlocu, zlocv, zlocw          ! temporary arrays
63      INTEGER  :: ji, jj, jk                                ! dummy loop indices
64
65     
66      IF( nn_diacfl == 1) THEN
67         IF( nn_timing == 1 )   CALL timing_start('dia_cfl')
68         ! setup timestep multiplier to account for initial Eulerian timestep
69         IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    dt = rdt
70         ELSE                                        ;    dt = rdt * 2.0
71         ENDIF
72
73             ! calculate Courant numbers
74         DO jk = 1, jpk
75            DO jj = 1, jpj
76               DO ji = 1, fs_jpim1   ! vector opt.
77
78                  ! Courant number for x-direction (zonal current)
79                  zcu_cfl(ji,jj,jk) = ABS(un(ji,jj,jk))*dt/e1u(ji,jj)
80
81                  ! Courant number for y-direction (meridional current)
82                  zcv_cfl(ji,jj,jk) = ABS(vn(ji,jj,jk))*dt/e2v(ji,jj)
83
84                  ! Courant number for z-direction (vertical current)
85                  zcw_cfl(ji,jj,jk) = ABS(wn(ji,jj,jk))*dt/fse3w_n(ji,jj,jk)
86               END DO
87            END DO         
88         END DO
89
90         ! calculate maximum values and locations
91         IF( lk_mpp ) THEN
92            CALL mpp_maxloc(zcu_cfl,umask,zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3))
93            CALL mpp_maxloc(zcv_cfl,vmask,zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3))
94            CALL mpp_maxloc(zcw_cfl,tmask,zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3))
95         ELSE
96            zlocu = MAXLOC( ABS( zcu_cfl(:,:,:) ) )
97            zcu_loc(1) = zlocu(1) + nimpp - 1
98            zcu_loc(2) = zlocu(2) + njmpp - 1
99            zcu_loc(3) = zlocu(3)
100            zcu_max = zcu_cfl(zcu_loc(1),zcu_loc(2),zcu_loc(3))
101
102            zlocv = MAXLOC( ABS( zcv_cfl(:,:,:) ) )
103            zcv_loc(1) = zlocv(1) + nimpp - 1
104            zcv_loc(2) = zlocv(2) + njmpp - 1
105            zcv_loc(3) = zlocv(3)
106            zcv_max = zcv_cfl(zcv_loc(1),zcv_loc(2),zcv_loc(3))
107
108            zlocw = MAXLOC( ABS( zcw_cfl(:,:,:) ) )
109            zcw_loc(1) = zlocw(1) + nimpp - 1
110            zcw_loc(2) = zlocw(2) + njmpp - 1
111            zcw_loc(3) = zlocw(3)
112            zcw_max = zcw_cfl(zcw_loc(1),zcw_loc(2),zcw_loc(3))
113         ENDIF
114     
115         ! write out to file
116         IF( lwp ) THEN
117            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)
118            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)
119            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)
120         ENDIF
121
122         ! update maximum Courant numbers from whole run if applicable
123         IF( zcu_max > cu_max ) THEN
124            cu_max = zcu_max
125            cu_loc = zcu_loc
126         ENDIF
127         IF( zcv_max > cv_max ) THEN
128            cv_max = zcv_max
129            cv_loc = zcv_loc
130         ENDIF
131         IF( zcw_max > cw_max ) THEN
132            cw_max = zcw_max
133            cw_loc = zcw_loc
134         ENDIF
135
136         ! at end of run output max Cu and Cv and close ascii file
137         IF( kt == nitend .AND. lwp ) THEN
138            ! to ascii file
139            WRITE(numcfl,*) '******************************************'
140            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)
141            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
142            WRITE(numcfl,*) '******************************************'
143            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)
144            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
145            WRITE(numcfl,*) '******************************************'
146            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)
147            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
148            CLOSE( numcfl ) 
149
150            ! to ocean output
151            WRITE(numout,*)
152            WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:'
153            WRITE(numout,*) '~~~~~~~~~~~~'
154            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) = (', 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) = (', cv_loc(1), cv_loc(2), cv_loc(3), ')'
157            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
158            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) = (', cw_loc(1), cw_loc(2), cw_loc(3), ')'
159            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
160
161         ENDIF
162
163         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl')
164      ENDIF
165
166   END SUBROUTINE dia_cfl
167
168   SUBROUTINE dia_cfl_init
169      !!----------------------------------------------------------------------
170      !!                  ***  ROUTINE dia_cfl_init  ***
171      !!                   
172      !! ** Purpose :   create output file, initialise arrays
173      !!----------------------------------------------------------------------
174
175
176      IF( nn_diacfl == 1 ) THEN
177         IF( nn_timing == 1 )   CALL timing_start('dia_cfl_init')
178
179         cu_max=0.0
180         cv_max=0.0
181         cw_max=0.0
182
183         ALLOCATE( zcu_cfl(jpi, jpj, jpk), zcv_cfl(jpi, jpj, jpk), zcw_cfl(jpi, jpj, jpk) )
184
185         zcu_cfl(:,:,:)=0.0
186         zcv_cfl(:,:,:)=0.0
187         zcw_cfl(:,:,:)=0.0
188
189         IF( lwp ) THEN
190            WRITE(numout,*)
191            WRITE(numout,*) 'dia_cfl     : Outputting CFL diagnostics to '//TRIM(clname)
192            WRITE(numout,*) '~~~~~~~~~~~~'
193            WRITE(numout,*)
194
195            ! create output ascii file
196            CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
197            WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k'
198            WRITE(numcfl,*) '******************************************'
199         ENDIF
200
201         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl_init')
202
203      ENDIF
204
205   END SUBROUTINE dia_cfl_init
206
207END MODULE diacfl
Note: See TracBrowser for help on using the repository browser.