[4715] | 1 | MODULE 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 | |
---|
| 45 | CONTAINS |
---|
| 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 | |
---|
| 210 | END MODULE diacfl |
---|