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_diagnostic/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2014/dev_r4650_UKMO12_CFL_diagnostic/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90 @ 4706

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

Added diacfl.F90 (copied from a vn3.4 local branch)
Modified step and step_oce to call it.
Currently uses a CPP key - this will be replaced by namelist control
File size: 8.7 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   !!
8   !!----------------------------------------------------------------------
9#if defined key_diacfl
10   !!----------------------------------------------------------------------
11   !!   dia_cfl        : Compute and output Courant numbers at each timestep
12   !!----------------------------------------------------------------------
13   USE oce             ! ocean dynamics and active tracers
14   USE dom_oce         ! ocean space and time domain
15   USE lib_mpp         ! distribued memory computing
16   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
17   USE in_out_manager  ! I/O manager
18   USE domvvl
19
20   IMPLICIT NONE
21   PRIVATE
22
23   LOGICAL, PUBLIC, PARAMETER   ::   lk_diacfl = .TRUE.    ! CFL diag flag
24   REAL(wp) :: cu_max, cv_max, cw_max                      ! daily max U Courant number
25   INTEGER, DIMENSION(3) :: cu_loc, cv_loc, cw_loc         ! daily max locations
26
27   INTEGER  :: numcfl                                       ! outfile unit
28   CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"      ! ascii filename
29
30   PUBLIC   dia_cfl       ! routine called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
37   !! $Id$
38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40
41
42CONTAINS
43
44
45   SUBROUTINE dia_cfl ( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE dia_cfl  ***
48      !!
49      !! ** Purpose :  Compute the Courant numbers Cu=u*dt/dx and Cv=v*dt/dy
50      !!               and output to ascii file 'cfl_diagnostics.ascii'
51      !!----------------------------------------------------------------------
52
53      INTEGER, INTENT(in) ::  kt                            ! ocean time-step index
54
55      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcu_cfl           ! Courant number arrays
56      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcv_cfl           !     "     "     "
57      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zcw_cfl           !     "     "     "
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( kt == nit000  ) THEN
67         cu_max=0.0
68         cv_max=0.0
69         cw_max=0.0
70         zcu_cfl(:,:,:)=0.0
71         zcv_cfl(:,:,:)=0.0
72         zcw_cfl(:,:,:)=0.0
73
74         IF( lwp ) THEN
75            WRITE(numout,*)
76            WRITE(numout,*) 'dia_cfl     : Outputting CFL diagnostics to '//TRIM(clname)
77            WRITE(numout,*) '~~~~~~~~~~~~'
78            WRITE(numout,*)
79
80            ! create output ascii file
81            CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
82            WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k'
83            WRITE(numcfl,*) '******************************************'
84         ENDIF
85      ENDIF
86
87      ! setup timestep multiplier to account for initial Eulerian timestep
88      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    dt = rdt
89      ELSE                                        ;    dt = rdt * 2.0
90      ENDIF
91
92      ! calculate Courant numbers
93      DO jk = 1, jpk
94         DO jj = 1, jpj
95            DO ji = 1, fs_jpim1   ! vector opt.
96
97               ! Courant number for x-direction (zonal current)
98               zcu_cfl(ji,jj,jk) = ABS(un(ji,jj,jk))*dt/e1u(ji,jj)
99
100               ! Courant number for y-direction (meridional current)
101               zcv_cfl(ji,jj,jk) = ABS(vn(ji,jj,jk))*dt/e2v(ji,jj)
102
103               ! Courant number for z-direction (vertical current)
104               zcw_cfl(ji,jj,jk) = ABS(wn(ji,jj,jk))*dt/fse3w_n(ji,jj,jk)
105            END DO
106         END DO         
107      END DO
108      CALL lbc_lnk( zcu_cfl(:,:,:), 'U', 1. )
109      CALL lbc_lnk( zcv_cfl(:,:,:), 'V', 1. )
110      CALL lbc_lnk( zcw_cfl(:,:,:), 'W', 1. )
111
112      ! calculate maximum values and locations
113      IF( lk_mpp ) THEN
114         CALL mpp_maxloc(zcu_cfl,umask,zcu_max, zcu_loc(1), zcu_loc(2), zcu_loc(3))
115         CALL mpp_maxloc(zcv_cfl,vmask,zcv_max, zcv_loc(1), zcv_loc(2), zcv_loc(3))
116         CALL mpp_maxloc(zcw_cfl,tmask,zcw_max, zcw_loc(1), zcw_loc(2), zcw_loc(3))
117      ELSE
118         zlocu = MAXLOC( ABS( zcu_cfl(:,:,:) ) )
119         zcu_loc(1) = zlocu(1) + nimpp - 1
120         zcu_loc(2) = zlocu(2) + njmpp - 1
121         zcu_loc(3) = zlocu(3)
122         zcu_max = zcu_cfl(zcu_loc(1),zcu_loc(2),zcu_loc(3))
123
124         zlocv = MAXLOC( ABS( zcv_cfl(:,:,:) ) )
125         zcv_loc(1) = zlocv(1) + nimpp - 1
126         zcv_loc(2) = zlocv(2) + njmpp - 1
127         zcv_loc(3) = zlocv(3)
128         zcv_max = zcv_cfl(zcv_loc(1),zcv_loc(2),zcv_loc(3))
129
130         zlocw = MAXLOC( ABS( zcw_cfl(:,:,:) ) )
131         zcw_loc(1) = zlocw(1) + nimpp - 1
132         zcw_loc(2) = zlocw(2) + njmpp - 1
133         zcw_loc(3) = zlocw(3)
134         zcw_max = zcw_cfl(zcw_loc(1),zcw_loc(2),zcw_loc(3))
135      ENDIF
136   
137      ! write out to file
138      IF( lwp ) THEN
139         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)
140         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)
141         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)
142      ENDIF
143
144      ! update maximum daily Courant numbers if applicable
145      IF( zcu_max > cu_max ) THEN
146         cu_max = zcu_max
147         cu_loc = zcu_loc
148      ENDIF
149      IF( zcv_max > cv_max ) THEN
150         cv_max = zcv_max
151         cv_loc = zcv_loc
152      ENDIF
153      IF( zcw_max > cw_max ) THEN
154         cw_max = zcw_max
155         cw_loc = zcw_loc
156      ENDIF
157
158      ! at end of run output max Cu and Cv and close ascii file
159      IF( kt == nitend .AND. lwp ) THEN
160         ! to ascii file
161         WRITE(numcfl,*) '******************************************'
162         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)
163         WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
164         WRITE(numcfl,*) '******************************************'
165         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)
166         WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
167         WRITE(numcfl,*) '******************************************'
168         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)
169         WRITE(numcfl,FMT='(3x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
170         CLOSE( numcfl ) 
171
172         ! to ocean output
173         WRITE(numout,*)
174         WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:'
175         WRITE(numout,*) '~~~~~~~~~~~~'
176         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), ')'
177         WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
178         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), ')'
179         WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
180         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), ')'
181         WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
182
183      ENDIF
184
185   END SUBROUTINE dia_cfl
186
187#else
188   !!----------------------------------------------------------------------
189   !!   Default option :                                       Dummy Module
190   !!----------------------------------------------------------------------
191   LOGICAL, PUBLIC, PARAMETER   ::   lk_diacfl = .FALSE.   ! CFL diag flag
192CONTAINS
193   SUBROUTINE dia_cfl( kt )        ! Empty routine
194      WRITE(*,*) 'dia_cfl: : You should not have seen this print! error?', kt
195   END SUBROUTINE dia_cfl
196#endif
197
198   !!======================================================================
199END MODULE diacfl
Note: See TracBrowser for help on using the repository browser.