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.
diacfl.F90 in branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

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

Last change on this file since 4715 was 4715, checked in by timgraham, 10 years ago

Added nn_diacfl to the namctl read in other configurations.
File size: 9.2 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         CALL lbc_lnk( zcu_cfl(:,:,:), 'U', 1. )
90         CALL lbc_lnk( zcv_cfl(:,:,:), 'V', 1. )
91         CALL lbc_lnk( zcw_cfl(:,:,:), 'W', 1. )
92
93         ! calculate maximum values and locations
94         IF( lk_mpp ) THEN
95            CALL mpp_maxloc(zcu_cfl,umask,zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3))
96            CALL mpp_maxloc(zcv_cfl,vmask,zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3))
97            CALL mpp_maxloc(zcw_cfl,tmask,zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3))
98         ELSE
99            zlocu = MAXLOC( ABS( zcu_cfl(:,:,:) ) )
100            zcu_loc(1) = zlocu(1) + nimpp - 1
101            zcu_loc(2) = zlocu(2) + njmpp - 1
102            zcu_loc(3) = zlocu(3)
103            zcu_max = zcu_cfl(zcu_loc(1),zcu_loc(2),zcu_loc(3))
104
105            zlocv = MAXLOC( ABS( zcv_cfl(:,:,:) ) )
106            zcv_loc(1) = zlocv(1) + nimpp - 1
107            zcv_loc(2) = zlocv(2) + njmpp - 1
108            zcv_loc(3) = zlocv(3)
109            zcv_max = zcv_cfl(zcv_loc(1),zcv_loc(2),zcv_loc(3))
110
111            zlocw = MAXLOC( ABS( zcw_cfl(:,:,:) ) )
112            zcw_loc(1) = zlocw(1) + nimpp - 1
113            zcw_loc(2) = zlocw(2) + njmpp - 1
114            zcw_loc(3) = zlocw(3)
115            zcw_max = zcw_cfl(zcw_loc(1),zcw_loc(2),zcw_loc(3))
116         ENDIF
117     
118         ! write out to file
119         IF( lwp ) THEN
120            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)
121            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)
122            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)
123         ENDIF
124
125         ! update maximum Courant numbers from whole run if applicable
126         IF( zcu_max > cu_max ) THEN
127            cu_max = zcu_max
128            cu_loc = zcu_loc
129         ENDIF
130         IF( zcv_max > cv_max ) THEN
131            cv_max = zcv_max
132            cv_loc = zcv_loc
133         ENDIF
134         IF( zcw_max > cw_max ) THEN
135            cw_max = zcw_max
136            cw_loc = zcw_loc
137         ENDIF
138
139         ! at end of run output max Cu and Cv and close ascii file
140         IF( kt == nitend .AND. lwp ) THEN
141            ! to ascii file
142            WRITE(numcfl,*) '******************************************'
143            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)
144            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
145            WRITE(numcfl,*) '******************************************'
146            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)
147            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
148            WRITE(numcfl,*) '******************************************'
149            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)
150            WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
151            CLOSE( numcfl ) 
152
153            ! to ocean output
154            WRITE(numout,*)
155            WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:'
156            WRITE(numout,*) '~~~~~~~~~~~~'
157            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), ')'
158            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
159            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), ')'
160            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
161            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), ')'
162            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
163
164         ENDIF
165
166         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl')
167      ENDIF
168
169   END SUBROUTINE dia_cfl
170
171   SUBROUTINE dia_cfl_init
172      !!----------------------------------------------------------------------
173      !!                  ***  ROUTINE dia_cfl_init  ***
174      !!                   
175      !! ** Purpose :   create output file, initialise arrays
176      !!----------------------------------------------------------------------
177
178
179      IF( nn_diacfl == 1 ) THEN
180         IF( nn_timing == 1 )   CALL timing_start('dia_cfl_init')
181
182         cu_max=0.0
183         cv_max=0.0
184         cw_max=0.0
185
186         ALLOCATE( zcu_cfl(jpi, jpj, jpk), zcv_cfl(jpi, jpj, jpk), zcw_cfl(jpi, jpj, jpk) )
187
188         zcu_cfl(:,:,:)=0.0
189         zcv_cfl(:,:,:)=0.0
190         zcw_cfl(:,:,:)=0.0
191
192         IF( lwp ) THEN
193            WRITE(numout,*)
194            WRITE(numout,*) 'dia_cfl     : Outputting CFL diagnostics to '//TRIM(clname)
195            WRITE(numout,*) '~~~~~~~~~~~~'
196            WRITE(numout,*)
197
198            ! create output ascii file
199            CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
200            WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k'
201            WRITE(numcfl,*) '******************************************'
202         ENDIF
203
204         IF( nn_timing == 1 )   CALL timing_stop('dia_cfl_init')
205
206      ENDIF
207
208   END SUBROUTINE dia_cfl_init
209
210END MODULE diacfl
Note: See TracBrowser for help on using the repository browser.