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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DIA/diacfl.F90 @ 11053

Last change on this file since 11053 was 10965, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : DIA and stpctl.F90. Just testing in ORCA1 so far.
  • Property svn:keywords set to Id
File size: 8.4 KB
Line 
1MODULE diacfl
2   !!======================================================================
3   !!                       ***  MODULE  diacfl  ***
4   !! Output CFL diagnostics to ascii file
5   !!======================================================================
6   !! History :  3.4  !  2010-03  (E. Blockley)  Original code
7   !!            3.6  !  2014-06  (T. Graham)  Removed CPP key & Updated to vn3.6
8   !!            4.0  !  2017-09  (G. Madec)  style + comments
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 domvvl          !
15   !
16   USE lib_mpp         ! distribued memory computing
17   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
18   USE in_out_manager  ! I/O manager
19   USE timing          ! Performance output
20
21   IMPLICIT NONE
22   PRIVATE
23
24   CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii"    ! ascii filename
25   INTEGER           :: numcfl                            ! outfile unit
26   !
27   INTEGER, DIMENSION(3) ::   nCu_loc, nCv_loc, nCw_loc   ! U, V, and W run max locations in the global domain
28   REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number
29
30!!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc !
31!!gm          I don't understand why.
32   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
33!!gm end
34
35   PUBLIC   dia_cfl       ! routine called by step.F90
36   PUBLIC   dia_cfl_init  ! routine called by nemogcm
37
38   !! * Substitutions
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
42   !! $Id$
43   !! Software governed by the CeCILL license (see ./LICENSE)
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE dia_cfl ( kt, Kmm )
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      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
55      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
56      !
57      INTEGER                ::   ji, jj, jk                            ! dummy loop indices
58      REAL(wp)               ::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars
59      INTEGER , DIMENSION(3) ::   iloc_u , iloc_v , iloc_w , iloc       ! workspace
60!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
61      !!----------------------------------------------------------------------
62      !
63      IF( ln_timing )   CALL timing_start('dia_cfl')
64      !
65      !                       ! setup timestep multiplier to account for initial Eulerian timestep
66      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2dt = rdt
67      ELSE                                        ;    z2dt = rdt * 2._wp
68      ENDIF
69      !
70      !               
71      DO jk = 1, jpk       ! calculate Courant numbers
72         DO jj = 1, jpj
73            DO ji = 1, fs_jpim1   ! vector opt.
74               zCu_cfl(ji,jj,jk) = ABS( uu(ji,jj,jk,Kmm) ) * z2dt / e1u  (ji,jj)      ! for i-direction
75               zCv_cfl(ji,jj,jk) = ABS( vv(ji,jj,jk,Kmm) ) * z2dt / e2v  (ji,jj)      ! for j-direction
76               zCw_cfl(ji,jj,jk) = ABS( ww(ji,jj,jk) ) * z2dt / e3w(ji,jj,jk,Kmm)   ! for k-direction
77            END DO
78         END DO         
79      END DO
80      !
81      !                    ! calculate maximum values and locations
82      IF( lk_mpp ) THEN
83         CALL mpp_maxloc( 'diacfl', zCu_cfl, umask, zCu_max, iloc_u )
84         CALL mpp_maxloc( 'diacfl', zCv_cfl, vmask, zCv_max, iloc_v )
85         CALL mpp_maxloc( 'diacfl', zCw_cfl, wmask, zCw_max, iloc_w )
86      ELSE
87         iloc = MAXLOC( ABS( zcu_cfl(:,:,:) ) )
88         iloc_u(1) = iloc(1) + nimpp - 1
89         iloc_u(2) = iloc(2) + njmpp - 1
90         iloc_u(3) = iloc(3)
91         zCu_max = zCu_cfl(iloc(1),iloc(2),iloc(3))
92         !
93         iloc = MAXLOC( ABS( zcv_cfl(:,:,:) ) )
94         iloc_v(1) = iloc(1) + nimpp - 1
95         iloc_v(2) = iloc(2) + njmpp - 1
96         iloc_v(3) = iloc(3)
97         zCv_max = zCv_cfl(iloc(1),iloc(2),iloc(3))
98         !
99         iloc = MAXLOC( ABS( zcw_cfl(:,:,:) ) )
100         iloc_w(1) = iloc(1) + nimpp - 1
101         iloc_w(2) = iloc(2) + njmpp - 1
102         iloc_w(3) = iloc(3)
103         zCw_max = zCw_cfl(iloc(1),iloc(2),iloc(3))
104      ENDIF
105      !
106      !                    ! write out to file
107      IF( lwp ) THEN
108         WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
109         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
110         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
111      ENDIF
112      !
113      !                    ! update maximum Courant numbers from whole run if applicable
114      IF( zCu_max > rCu_max ) THEN   ;   rCu_max = zCu_max   ;   nCu_loc(:) = iloc_u(:)   ;   ENDIF
115      IF( zCv_max > rCv_max ) THEN   ;   rCv_max = zCv_max   ;   nCv_loc(:) = iloc_v(:)   ;   ENDIF
116      IF( zCw_max > rCw_max ) THEN   ;   rCw_max = zCw_max   ;   nCw_loc(:) = iloc_w(:)   ;   ENDIF
117
118      !                    ! at end of run output max Cu and Cv and close ascii file
119      IF( kt == nitend .AND. lwp ) THEN
120         ! to ascii file
121         WRITE(numcfl,*) '******************************************'
122         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3)
123         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCu_max
124         WRITE(numcfl,*) '******************************************'
125         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3)
126         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCv_max
127         WRITE(numcfl,*) '******************************************'
128         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3)
129         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCw_max
130         CLOSE( numcfl ) 
131         !
132         ! to ocean output
133         WRITE(numout,*)
134         WRITE(numout,*) 'dia_cfl : Maximum Courant number information for the run '
135         WRITE(numout,*) '~~~~~~~'
136         WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', z2dt/rCu_max
137         WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', z2dt/rCv_max
138         WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', z2dt/rCw_max
139      ENDIF
140      !
141      IF( ln_timing )   CALL timing_stop('dia_cfl')
142      !
143   END SUBROUTINE dia_cfl
144
145
146   SUBROUTINE dia_cfl_init
147      !!----------------------------------------------------------------------
148      !!                  ***  ROUTINE dia_cfl_init  ***
149      !!                   
150      !! ** Purpose :   create output file, initialise arrays
151      !!----------------------------------------------------------------------
152      !
153      IF(lwp) THEN
154         WRITE(numout,*)
155         WRITE(numout,*) 'dia_cfl : Outputting CFL diagnostics to ',TRIM(clname), ' file'
156         WRITE(numout,*) '~~~~~~~'
157         WRITE(numout,*)
158         !
159         ! create output ascii file
160         CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
161         WRITE(numcfl,*) 'Timestep  Direction  Max C     i    j    k'
162         WRITE(numcfl,*) '******************************************'
163      ENDIF
164      !
165      rCu_max = 0._wp
166      rCv_max = 0._wp
167      rCw_max = 0._wp
168      !
169!!gm required to work
170      ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) )
171!!gm end
172      !     
173   END SUBROUTINE dia_cfl_init
174
175   !!======================================================================
176END MODULE diacfl
Note: See TracBrowser for help on using the repository browser.