- Timestamp:
- 2017-12-01T18:44:09+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90
r8329 r8882 1 1 MODULE diacfl 2 !!====================================================================== ========2 !!====================================================================== 3 3 !! *** MODULE diacfl *** 4 4 !! Output CFL diagnostics to ascii file 5 !!====================================================================== ========6 !! History : 1.0! 2010-03 (E. Blockley) Original code7 !! ! 2014-06 (T Graham)Removed CPP key & Updated to vn3.68 !! 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 9 !!---------------------------------------------------------------------- 10 10 !! dia_cfl : Compute and output Courant numbers at each timestep … … 12 12 USE oce ! ocean dynamics and active tracers 13 13 USE dom_oce ! ocean space and time domain 14 USE domvvl ! 15 ! 14 16 USE lib_mpp ! distribued memory computing 15 17 USE lbclnk ! ocean lateral boundary condition (or mpp link) 16 18 USE in_out_manager ! I/O manager 17 USE domvvl18 19 USE timing ! Performance output 19 20 … … 21 22 PRIVATE 22 23 23 REAL(wp) :: cu_max, cv_max, cw_max ! Run max U Courant number24 INTEGER , DIMENSION(3) :: cu_loc, cv_loc, cw_loc ! Run max locations25 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcu_cfl ! Courant number arrays26 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcv_cfl ! Courant number arrays27 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zcw_cfl ! Courant number arrays24 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 28 29 29 INTEGER :: numcfl ! outfile unit 30 CHARACTER(LEN=50) :: clname="cfl_diagnostics.ascii" ! ascii filename 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 31 34 32 35 PUBLIC dia_cfl ! routine called by step.F90 … … 40 43 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 41 44 !!---------------------------------------------------------------------- 42 43 44 45 CONTAINS 45 46 46 47 47 SUBROUTINE dia_cfl ( kt ) … … 52 52 !! and output to ascii file 'cfl_diagnostics.ascii' 53 53 !!---------------------------------------------------------------------- 54 INTEGER, INTENT(in) :: kt ! ocean time-step index 55 ! 56 INTEGER :: ji, jj, jk ! dummy loop indices 57 REAL(wp):: z2dt, zCu_max, zCv_max, zCw_max ! local scalars 58 INTEGER , DIMENSION(3) :: iloc_u , iloc_v , iloc_w , iloc ! workspace 59 !!gm this does not work REAL(wp), DIMENSION(jpi,jpj,jpk) :: zCu_cfl, zCv_cfl, zCw_cfl ! workspace 60 !!---------------------------------------------------------------------- 61 ! 62 IF( nn_timing == 1 ) CALL timing_start('dia_cfl') 63 ! 64 ! ! setup timestep multiplier to account for initial Eulerian timestep 65 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2dt = rdt 66 ELSE ; z2dt = rdt * 2._wp 67 ENDIF 68 ! 69 ! 70 DO jk = 1, jpk ! calculate Courant numbers 71 DO jj = 1, jpj 72 DO ji = 1, fs_jpim1 ! vector opt. 73 zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u (ji,jj) ! for i-direction 74 zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v (ji,jj) ! for j-direction 75 zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk) ! for k-direction 76 END DO 77 END DO 78 END DO 79 ! 80 ! ! calculate maximum values and locations 81 IF( lk_mpp ) THEN 82 CALL mpp_maxloc( zCu_cfl, umask, zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) ) 83 CALL mpp_maxloc( zCv_cfl, vmask, zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) ) 84 CALL mpp_maxloc( zCw_cfl, wmask, zCw_max, iloc_w(1), iloc_w(2), iloc_w(3) ) 85 ELSE 86 iloc = MAXLOC( ABS( zcu_cfl(:,:,:) ) ) 87 iloc_u(1) = iloc(1) + nimpp - 1 88 iloc_u(2) = iloc(2) + njmpp - 1 89 iloc_u(3) = iloc(3) 90 zCu_max = zCu_cfl(iloc(1),iloc(2),iloc(3)) 91 ! 92 iloc = MAXLOC( ABS( zcv_cfl(:,:,:) ) ) 93 iloc_v(1) = iloc(1) + nimpp - 1 94 iloc_v(2) = iloc(2) + njmpp - 1 95 iloc_v(3) = iloc(3) 96 zCv_max = zCv_cfl(iloc(1),iloc(2),iloc(3)) 97 ! 98 iloc = MAXLOC( ABS( zcw_cfl(:,:,:) ) ) 99 iloc_w(1) = iloc(1) + nimpp - 1 100 iloc_w(2) = iloc(2) + njmpp - 1 101 iloc_w(3) = iloc(3) 102 zCw_max = zCw_cfl(iloc(1),iloc(2),iloc(3)) 103 ENDIF 104 ! 105 ! ! write out to file 106 IF( lwp ) THEN 107 WRITE(numcfl,FMT='(2x,i4,5x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3) 108 WRITE(numcfl,FMT='(11x, a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3) 109 WRITE(numcfl,FMT='(11x, a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') 'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3) 110 ENDIF 111 ! 112 ! ! update maximum Courant numbers from whole run if applicable 113 IF( zCu_max > rCu_max ) THEN ; rCu_max = zCu_max ; nCu_loc(:) = iloc_u(:) ; ENDIF 114 IF( zCv_max > rCv_max ) THEN ; rCv_max = zCv_max ; nCv_loc(:) = iloc_v(:) ; ENDIF 115 IF( zCw_max > rCw_max ) THEN ; rCw_max = zCw_max ; nCw_loc(:) = iloc_w(:) ; ENDIF 54 116 55 INTEGER, INTENT(in) :: kt ! ocean time-step index 117 ! ! at end of run output max Cu and Cv and close ascii file 118 IF( kt == nitend .AND. lwp ) THEN 119 ! to ascii file 120 WRITE(numcfl,*) '******************************************' 121 WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3) 122 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCu_max 123 WRITE(numcfl,*) '******************************************' 124 WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3) 125 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCv_max 126 WRITE(numcfl,*) '******************************************' 127 WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3) 128 WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCw_max 129 CLOSE( numcfl ) 130 ! 131 ! to ocean output 132 WRITE(numout,*) 133 WRITE(numout,*) 'dia_cfl : Maximum Courant number information for the run ' 134 WRITE(numout,*) '~~~~~~~' 135 WRITE(numout,*) ' Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', z2dt/rCu_max 136 WRITE(numout,*) ' Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', z2dt/rCv_max 137 WRITE(numout,*) ' Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', z2dt/rCw_max 138 ENDIF 139 ! 140 IF( nn_timing == 1 ) CALL timing_stop('dia_cfl') 141 ! 142 END SUBROUTINE dia_cfl 56 143 57 REAL(wp) :: zcu_max, zcv_max, zcw_max ! max Courant numbers per timestep58 INTEGER, DIMENSION(3) :: zcu_loc, zcv_loc, zcw_loc ! max Courant number locations59 60 REAL(wp) :: dt ! temporary scalars61 INTEGER, DIMENSION(3) :: zlocu, zlocv, zlocw ! temporary arrays62 INTEGER :: ji, jj, jk ! dummy loop indices63 64 65 IF( nn_diacfl == 1) THEN66 IF( nn_timing == 1 ) CALL timing_start('dia_cfl')67 ! setup timestep multiplier to account for initial Eulerian timestep68 IF( neuler == 0 .AND. kt == nit000 ) THEN ; dt = rdt69 ELSE ; dt = rdt * 2.070 ENDIF71 72 ! calculate Courant numbers73 DO jk = 1, jpk74 DO jj = 1, jpj75 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 DO86 END DO87 END DO88 89 ! calculate maximum values and locations90 IF( lk_mpp ) THEN91 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 ELSE95 zlocu = MAXLOC( ABS( zcu_cfl(:,:,:) ) )96 zcu_loc(1) = zlocu(1) + nimpp - 197 zcu_loc(2) = zlocu(2) + njmpp - 198 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 - 1103 zcv_loc(2) = zlocv(2) + njmpp - 1104 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 - 1109 zcw_loc(2) = zlocw(2) + njmpp - 1110 zcw_loc(3) = zlocw(3)111 zcw_max = zcw_cfl(zcw_loc(1),zcw_loc(2),zcw_loc(3))112 ENDIF113 114 ! write out to file115 IF( lwp ) THEN116 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 ENDIF120 121 ! update maximum Courant numbers from whole run if applicable122 IF( zcu_max > cu_max ) THEN123 cu_max = zcu_max124 cu_loc = zcu_loc125 ENDIF126 IF( zcv_max > cv_max ) THEN127 cv_max = zcv_max128 cv_loc = zcv_loc129 ENDIF130 IF( zcw_max > cw_max ) THEN131 cw_max = zcw_max132 cw_loc = zcw_loc133 ENDIF134 135 ! at end of run output max Cu and Cv and close ascii file136 IF( kt == nitend .AND. lwp ) THEN137 ! to ascii file138 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 output150 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 ENDIF164 165 IF( nn_timing == 1 ) CALL timing_stop('dia_cfl')166 ENDIF167 168 END SUBROUTINE dia_cfl169 144 170 145 SUBROUTINE dia_cfl_init … … 174 149 !! ** Purpose : create output file, initialise arrays 175 150 !!---------------------------------------------------------------------- 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 151 ! 152 IF(lwp) THEN 153 WRITE(numout,*) 154 WRITE(numout,*) 'dia_cfl : Outputting CFL diagnostics to ',TRIM(clname), ' file' 155 WRITE(numout,*) '~~~~~~~' 156 WRITE(numout,*) 157 ! 158 ! create output ascii file 159 CALL ctl_opn( numcfl, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 ) 160 WRITE(numcfl,*) 'Timestep Direction Max C i j k' 161 WRITE(numcfl,*) '******************************************' 205 162 ENDIF 206 163 ! 164 rCu_max = 0._wp 165 rCv_max = 0._wp 166 rCw_max = 0._wp 167 ! 168 !!gm required to work 169 ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) ) 170 !!gm end 171 ! 207 172 END SUBROUTINE dia_cfl_init 208 173 174 !!====================================================================== 209 175 END MODULE diacfl
Note: See TracChangeset
for help on using the changeset viewer.
