program BARO c c This program performs the nondivegent barotropic c model forecast. The prediction equation uses the c the streamfunction as its predicted variable. c The Eulerian frame is used. the advective term c is computed through the 9 point Arakawa jacobian. c Time integration is performed through the Matsuno c time scheme. Some properties of the nondivergent c barotropic fluid for example the maximum absolute c vorticity and the domain mean kinetic energy c are also computed. c c Originator Florida State University,last revision c L. Bounoua ,1994 c c Input : c streamfunction on unite 23. this file is c read in subroutine INIT . c c Output c : (i) forecast streamfunctions at preset c intervals of nint hours (see comments on c symbols below) .These are written onto unit c 17 using . the first record are the c streamfunction fields at 0 hour . subsequent c records contain the forecast records c which are sequentially written according c to the forecast hour. c c (ii) comment statements on energy computations c and vorticity parameters during time integrati c on.These are written onto unit 6 and normally c this file is assigned to the screen. c c Definition : c c l : number of grid points in the zonal c direction. c m : number of grid points in the meridional c direction. c n : number of vertical levels (1 in this c case). c l1 : l-1 c m1 : m-1 c n1 : n-1 c l2 : l-2 c m2 : m-2 c ladd : number of grid points added to make c domain cyclic. c phi(m) : latitudes of each row of grid points c cor(m) : coriolis parameter for each row of grid c points. c dy : grid size in the meridional direction c in meters. c dx(m) : grid size in the zonal direction in c meters for each row of grid points. c dt : time step of integration in seconds. c nint : time interval in hours to output c forecast products. this is later c modified to represent the number c of time steps before the next c forecast products are output to file c 17. c tend : Total length in hours of the forecast. c asml(n) : smallest value of the absolute vorticity c alrg(n) : largest value of the absolute vorticity. c psi(l,m,n) : streamfunction field. c a1(l,m,n) : array containing either the relative c or absolute vorticity. c a2(l,m,n) : another working array containing either c the relative or absolute vorticity. c b1(l,m,n) : array containing the advective term, c jacobian of(psi,absolute vorticity). c c Subroutines c called : bound,energy,equal, init,jac, c lap,large,relaxt, small, vort c parameter( l = 51,m = 21,n = 1) c c common l1, m1, l2, m2, & ladd, phi(m), cor(m), dy, & dx(m), dt, nint, mint, & tend, time, asml(n), alrg(n), & psi(l,m,n),a1(l,m,n),a2(l,m,n),b1(l,m,n), & slat ,dphi c write(6,1000) c c Call subroutine INIT to read in initial parameters c and to define the initial state . c call INIT c c Write the initial state to output file . c do 10200 k = 1, n write (17,878) ((psi(i,j,k),i=1,l),j=1,m) 10200 continue 878 format(6e13.6) c c Counts the number of timesteps in nrst . c nrst = 0 c c The largest and smallest value of the absolute vorticity in c each time step are calculated by subroutines LARGE and c SMALL respectively . c call LARGE (a1,l,m,n,alrg) call SMALL (a1,l,m,n,asml) write(6,1001)nrst+1 write(6,1002) time,alrg(1),asml(1) 1001 format(10x,'number of time steps = ',i3) 1002 format(/,5x,'time =',f8.0,4x,'max,min abs vort =',e11.4, & 5x,e11.4,/) c c The energy parameters between 2 latitudes at time = 0 c is calculated using subroutine ENERGY. c jslat1 = 2 jnlat1 = m1 jslat2 = 6 jnlat2 = m-jslat2+1 c call ENERGY (psi,l,m,n,l1,dy,dx,time,jslat1,jnlat1,slat,dphi) call ENERGY (psi,l,m,n,l1,dy,dx,time,jslat2,jnlat2,slat,dphi) c c Time integration following the Matsuno time schem c is carried out. The scheme consists of two steps, c namely the predictor step and the corrector step. c c Predicteur Step : c c First the advective term i.e. the jacobian c of psi and delsquared of (psi + f ) is computed . c 13 call JACMOD (b1,psi,a1,dx,dy,l,m,n,l1,m1,l2,m2,ladd) c c The time integration for the predictor step is c carried out next. This gives the first guess relative c vorticity field. It is stored in array a2. c do 10202 j = 1, m do 10202 i = 1, l 10202 a2(i,j,1) = (a1(i,j,1)-b1(i,j,1)*dt-cor(j)) do 10204 i = 1, l a2(i,1,1) = 0. 10204 a2(i,m,1) = 0. c c Corrector Step : c c Relaxation of a2 is done to obtain the first guess psi field. c with this, the first guess absolute vorticity field (a2+f) , c and the first guess advective term,j(psi,delsquared (psi + f) c are computed . c call RELAXT (psi,a2,dx,dy,l,l,m,l1,m1,m2,1) call VORT (a2,cor,l,m,n) call JACMOD (b1,psi,a2,dx,dy,l,m,n,l1,m1,l2,m2,ladd) c c The time integration for the corrector step is carried out c here. the second guess relative vorticity field is stored c in array a2. c do 10206 j = 1, m do 10206 i = 1, l 10206 a2(i,j,1) = (a1(i,j,1)-b1(i,j,1)*dt-cor(j)) do 10208 i = 1, l a2(i,1,1) = 0. 10208 a2(i,m,1) = 0. c c a2 is relaxed to obtain the second guess psi field. the second c guess absolute vorticity field (a2+f), is calculated and c transfered into array a1 by subroutine equal for the next c time step computations. c call RELAXT (psi,a2,dx,dy,l,l,m,l1,m1,m2,1) call VORT (a2,cor,l,m,n) call EQUAL (a1,a2,l,m,n) c c The time integration counter, time, is increased by c timestep dt . c nrst = nrst + 1 time = time + dt write(6,1001)nrst+1 c c The largest and smallest value of the absolute vorticity in c each time step are calculated by subroutines LARGE and c SMALL respectively. c call LARGE (a2,l,m,n,alrg) call SMALL (a2,l,m,n,asml) write(6,1002) time,alrg(1),asml(1) c c The energy parameters between 2 latitudes for each time step c are calculated using subroutine ENERGY c call ENERGY (psi,l,m,n,l1,dy,dx,time,jslat1,jnlat1,slat,dphi) call ENERGY (psi,l,m,n,l1,dy,dx,time,jslat2,jnlat2,slat,dphi) c c The process leading to another time integration is c is repeated. However, if time = nint (the preset c output interval), the psi field is written into c the output file . If time = tend (the preset time c to end forecast),the program execution is terminated. c if (nrst-nint) 13,14,13 14 continue do 10210 k = 1, n 10210 write (17,878) ((psi(i,j,k),i=1,l),j=1,m) nint = nint+mint c print *,'time ,tend',time,tend if (time-tend) 13,16,16 16 continue 9998 continue 1000 format(/) stop end subroutine LARGE (a,l,m,n,al) c c This subroutine searches for the largest value c of the variable at each level of a 3 dimensional c field. the search is carried out over the domain c one grid point inwards from the southern and c northern boundaries. c c Definitions : c c a(l,m,n) : array containing variable whose largest c value at each level is to be determined. c l : number of grid points in the zonal c direction. c m : number of grid points in the meridional c direction. c n : number of levels in the vertical. c al(n) : largest value of variable at level k. c dimension a(l,m,n),al(n) c m1 = m-1 do 10250 k = 1,n al(k) = a(1,2,k) do 10250 j = 2,m1 do 10250 i = 1,l if (al(k)-a(i,j,k)) 10251,10250,10250 10251 al(k) = a(i,j,k) 10250 continue return end subroutine INIT c c The subroutine defines the initial parameters required by the c main program BARO .These include the time step of integration, c the number of hours of forecasts to be done, the forecast output c interval, the grid sizes and the initial streamfunction field. c c Definitions : c the symbols used in the common block c are similar to those in program BARO. c c slat : southern most latitude of integration c domain. c dphi : grid size in degrees latitude/longitude. c parameter(l=51,m=21,n=1) c common l1, m1, l2, m2, 1 ladd, phi(m), cor(m), dy, 2 dx(m), dt, nint, mint, 3 tend, time, asml(n), alrg(n), 4 psi(l,m,n),a1(l,m,n),a2(l,m,n), b1(l,m,n), 5 slat ,dphi c c Define control parameters c slat = 0.0 dphi = 2.5 dt = 2700. nint = 12 time = 0. tend = 24. ladd = 6 tend = tend*3600. nint = int(nint*3600./dt) mint = nint l1 = l-1 l2 = l-2 m1 = m-1 m2 = m-2 n1 = n-1 pi = 3.1415927 rad = pi/180 phi(1) = slat cor(1) = 2.*7.292e-5*sin(phi(1)*rad) do 10240 j = 2, m phi(j) = phi(j-1)+dphi cor(j) = 2.*7.292e-5*sin(phi(j)*rad) 10240 continue dy = dphi*111.1*1000. do 10241 j = 1, m dx(j) = dy*cos(phi(j)*rad) 10241 continue c open(23,file='infield_out.dat',status='old ') open(17,file='baro_out.dat ',status='unknown') c c Read initial field from unit 23 and compute c Laplacian and Vorticity of initial field. c read(23,879) ((psi(i,j,1),i=1,l),j=1,m) 879 format(6e13.6) c call LAPMOD (a1,psi,dx,dy,l,m,n,l1,m1,l2,m2,ladd) call VORT (a1,cor,l,m,n) return end subroutine JACMOD (a,b,c,dx,dy,l,m,n,l1,m1,l2,m2,ladd) c c This subroutine computes the advective term, c the Jacobian of (psi, delsquared psi + f), c following Arakawa (1966). a nine-point stencil is c used . Note that this subroutine is similar in pri- c nciples to subroutine JAC ,but the cyclicity condi- c tion has been added . c c Defiitions : c c a(l,m,n) : this array contains the Jacobian of c psi and delsquared psi + f, where f is c the earth's vorticity. c c b(l,m,n) : the streamfunction field,psi. c c c(l,m,n) : the absolute vorticity or delsquared c (psi + f). c dx(m) : grid size in the zonal direction for c each row of grid points in meters. c dy : grid size in the meridional direction c in meters. c l : number of grid points in the zonal c direction. c m : number of grid points in the c meridional direction. c n : number of levels in the vertical. c l1 : l-1 c m1 : m-1 c l2 : l-2 c m2 : m-2 c ladd : number of grid points added to the c eastern boundary to make the domain c cyclic in the zonal direction. c real a(l,m,n), b(l,m,n), c(l,m,n), dx(m) k = 1 do 10140 j = 2, m1 dm = 12.*dx(j)*dy jp1 = j+1 jm1 = j-1 do 10140 i = 1, l if (i-1)80,80,81 80 im1 = l1 ip1 = 2 go to 83 81 if(i-l)82,80,80 82 im1 = i-1 ip1 = i+1 83 continue a(i,j,k) =(b(i,j-1,k)+b(ip1,j-1,k)-b(i,j+1,k)-b(ip1,j+1,k)) & *(c(ip1,j,k)-c(i,j,k))+(b(im1,j-1,k)+b(i,j-1,k)-b(im1,j+1,k & )-b(i,j+1,k))*(c(i,j,k)-c(im1,j,k))+(b(ip1,j,k)+b(ip1,j+1,k & )-b(im1,j,k)-b(im1,j+1,k))*(c(i,j+1,k)-c(i,j,k))+(b(ip1,j-1 & ,k)+b(ip1,j,k)-b(im1,j-1,k)-b(im1,j,k))*(c(i,j,k)-c(i,j-1,k & ))+(b(ip1,j,k)-b(i,j+1,k))*(c(ip1,j+1,k)-c(i,j,k))+(b(i,j-1 & ,k)-b(im1,j,k))*(c(i,j,k)-c(im1,j-1,k))+(b(i,j+1,k)-b(im1,j & ,k))*(c(im1,j+1,k)-c(i,j,k))+(b(ip1,j,k)-b(i,j-1,k))*(c(i,j & ,k)-c(ip1,j-1,k)) c a(i,j,k) = a(i,j,k)/dm 10140 continue do 10141 i= 1, l if (i-1)70,70,71 70 im1 = l1 ip1 = 2 go to 73 71 if (i-l) 72,70,70 72 im1 = i-1 ip1 = i+1 73 continue c a(i,1,k) = (b(i,1,k)+b(ip1,1,k)-b(i,2,k)-b(ip1,2,k))*(c(i,1,k)+ & c(ip1,1,k))-(b(im1,1,k)+b(i,1,k)-b(im1,2,k)-b(i,2,k))*(c(im1 & ,1,k)+c(i,1,k))+(b(ip1,1,k)+b(ip1,2,k)-b(im1,1,k)-b(im1,2,k) & )*(c(i,1,k)+c(i,2,k))+(b(ip1,1,k)-b(i,2,k))*(c(i,1,k)+c(ip1, & 2,k))+(b(i,2,k)-b(im1,1,k))*(c(im1,2,k)+c(i,1,k)) c dm = 12.*dx(1)*dy a(i,1,k) = a(i,1,k)/dm a(i,m,k) = (b(i,m-1,k)+b(ip1,m-1,k)-b(i,m,k)-b(ip1,m,k))*(c(i,m, & k)+c(ip1,m,k))-(b(im1,m-1,k)+b(i,m-1,k)-b(im1,m,k)-b(i,m,k))* & (c(im1,m,k)+c(i,m,k))-(b(ip1,m-1,k)+b(ip1,m,k)-b(im1,m-1,k)-b & (im1,m,k))*(c(i,m-1,k)+c(i,m,k))-(b(i,m-1,k)-b(im1,m,k))*(c(i & m1,m-1,k)+c(i,m,k))-(b(ip1,m,k)-b(i,m-1,k))*(c(i,m,k)+c(ip1,m & -1,k)) c dm = 12.*dx(m)*dy a(i,m,k) = a(i,m,k)/dm 10141 continue c c The next do loop is for the non-cyclic case only. c if (ladd.eq.0) then do 10142 j = 1, m a(1,j,k) = a(2,j,k) a(l,j,k) = a(l1,j,k) 10142 continue endif return end subroutine EQUAL (a,b,l,m,n) c c This subroutine sets the values of the elements of c an array (a) equal to that of another array (b). c c Definition : c c a(l,m,n) : array in which the values of the c elements are to be set. c b(l,m,n) : array containing the values to set array c a(l,m,n). c l : number of grid points in the zonal c direction. c m : number of grid points in the meridional c direction. c n : number of levels in the vertical. c dimension a(l,m,n), b(l,m,n) k = 1 do 10230 j = 1, m do 10230 i = 1, l 10230 a(i,j,k) = b(i,j,k) return end subroutine VORT (a,cor,l,m,n) c c This subroutine computes the absolute vorticity c by adding the relative vorticity to the Coriolis c parameter. c c Definitions : c c a(l,m,n) : upon input, the array contains the relative c vorticity. Upon returns it contains the abs- c olute vorticity. c cor(m) : Coriolis parameter at each row of grid c points. c l : number of grid points in the zonal c direction. c m : number of grid points in the meridional c direction. c n : number of levels in the vertical. c real a(l,m,n), cor(m) k = 1 do 10280 j = 1, m do 10280 i = 1, l 10280 a(i,j,k) = a(i,j,k)+cor(j) return end subroutine ENERGY (psi,l,m,n,l1,dy,dx,time,j1,j2,slat,dphi) c c This subroutine calculates the domain average c (1) kinetic energy, (2) zonal kinetic energy, (3) c eddy kinetic energy and (4) the zonal to eddy c kinetic energy exchange. the domain is confined c to a latitudinal belt between grid rows j1 and j2. c See chapter 13 for details. c c Definition : c c psi(l,m,n) : streamfunction field. c l : number of grid points in the zonal c direction. c m : number of grid points in the meridional c direction. c n : number of levels in the vertical. c l1 : l-1 c dy : grid size in meters in the meridional c direction. c dx(m) : grid size in meters in the zonal c direction for each row of grid points. c time : time of forecast in seconds. c j1 : grid row number of southern boundary c of domain used for computing the energy c parameters. c j2 : grid row number of northern boundary of c domain used for computing the energy c parameters. c real psi(l,m,n),dx(m) c c Initialize parameters c ak = 0. akbar = 0. akzke = 0. c do 10221 j = j1, j2 jp1 = j+1 jm1 = j-1 ubar = 0 duvdy = 0 do 10220 i = 1,l1 im1 = i-1 if (i.eq.1) im1 = l1 u = -(psi(i,jp1,1)-psi(i,jm1,1))/(2.*dy) v = (psi(i+1,j,1)-psi(im1,j,1))/(2.*dx(j)) ak = ak+(u**2+v**2)/(2*l1*(j2-j1+1)) ubar = ubar+u/l1 uvn = -(psi(i,jp1,1)-psi(i,j,1))/dy*(psi(i+1,jp1,1) & -psi(im1,jp1,1)+psi(i+1,j,1)-psi(im1,j,1)) & /(4.*dx(j)) uvs = -(psi(i,j,1)-psi(i,jm1,1))/dy*(psi(i+1,j,1) & -psi(im1,j,1)+psi(i+1,jm1,1)-psi(im1,jm1,1)) & /(4.*dx(j)) duvdy = (uvn-uvs)/(dy*l1)+duvdy 10220 continue akbar = akbar+ubar**2/2./(j2-j1+1) akzke = akzke+ubar*duvdy/(j2-j1+1) 10221 continue akprim = ak-akbar alats = slat + (j1-1)*dphi alatn = slat + (j2-1)*dphi write (6,1000)alats,alatn 1000 format(5x,'energy exchange between ,lat=',f5.2,' and lat=',f5.2) write(6,1001) ak,akbar,akprim ,akzke 1001 format(5x,'k ,kbar, kprime, kzke ',4(1x,e11.4),/) return end subroutine SMALL (a,l,m,n,al) c c This subroutine searches for the smallest value c of the variable at each level of a 3 dimensional c field. the search is carried out over the domain c one grid point inwards from the southern and c northern boundaries. c c Definitions : c c a(l,m,n) : array containing variable whose smallest c value at each level is to be determined. c l : number of grid points in the zonal c direction. c m : number of grid points in the meridional c direction. c n : number of levels in the vertical. c al(n) : smallest value of variable at level k. c c real a(l,m,n), al(n) m1 = m-1 do 10270 k = 1,n al(k) = a(1,2,k) do 10270 j = 2,m1 do 10270 i = 1,l if (a(i,j,k)-al(k)) 10271,10270,10270 10271 al(k) = a(i,j,k) 10270 continue return end subroutine RELAXT (x,y,dx,dy,ll,l,m,l1,m1,m2,iskip) c c This subroutine solves the poisson equation, c delsquared (x) = y, through the sequential c over-relaxation method . This method assumes that c the boundary values of x are known. in particular, c this subroutine would be used to solve the equation, c delsquared ( psi) = zeta 'relative vorticity'. c c Definitions : c c x(l,m,n) : streamfunction field. c y(l,m,n) : relative vorticity. c dx(m) : grid size in meters in th e zonal c direction. c dy : grid size in meters in the meridional c direction . c ll : number of grid points in the zonal c direction. c l : (l-2) defines the number of grid c points in the zonal direction to be c relaxed. c m : number of grid points in the meridional c direction. c l1 : l-1 c m1 : m-1 c m2 : m-2 c iskip : control integer which determines the c values the first guess field of x(l,m,n) c assumes . if (iskip .le. 0) , a zero c guess field is assumed, otherwise it c takes the previous values of x(l,m,n) c as its guess field. c c eps : the maximum tolerance error for psi c real x(ll,m),y(ll,m),dx(m) write(6,1000) 1000 format(//) eps = 1000. c c 'nsc' is the counter for number of scans (iterations) done. c 'ia' is the maximum number of scans allowed and 'alfa' is c the (relaxation coefficient)/4. nsc = 0 ia = 500 lsc = -1 alfa = 0.45 if (iskip.le.0) then do 10260 j = 2, m1 do 10260 i = 2, l1 x(i,j) = 0. 10260 continue endif npts = (l-2)*(m-2) 15 nrel = 0 do 10261 j = 2, m1 jp1 = j+1 jm1 = j-1 do 10261 i = 2,l1 if (i-l) 5,6,6 5 if (i-1) 6,6,7 6 r = (x(2,j)+x(l1,j)-2.*x(i,j))/dx(j)**2 & +(x(i,jp1)+x(i,jm1)-2.*x(i,j))/dy**2 r = (r-y(i,j))*dy*dx(j) go to 21 7 r = (x(i+1,j)+x(i-1,j)-2.*x(i,j))/dx(j)**2 & +(x(i,jp1)+x(i,jm1)-2.*x(i,j))/dy**2 r = (r-y(i,j))*dy*dx(j) 21 if (lsc-nsc) 29,29,30 29 x(i,j)=x(i,j)+alfa*r 30 if (abs(r).le.eps) nrel = nrel+1 10261 continue nsc = nsc + 1 if (nrel-npts) 13,14,14 14 if (lsc .ge. nsc) go to 300 18 lsc = nsc + 1 13 if (nsc.lt.ia) go to 15 300 continue return end subroutine LAPMOD (a,b,dx,dy,l,m,n,l1,m1,l2,m2,ladd) c c This subroutine LAP calculates the Laplacian of any c scalar field. the standard second order finite c differencing is used . two forms of boundary c conditions are included, one that is cyclic in c the zonal direction, and the other that is non c cyclic in the zonal direction. for the cyclic c case, the southern and northern boundary values are c formed by linear extrapolation outwards , while for c the non-cyclic case , the boundary values are formed c by extending the values one grid point outwards. c Note that this subroutine is similar in principles c to subroutine JAC ,but the cyclicity condition has be- c en added . c c Definitions : c c a(l,m,n) : Laplacian of b. for our purpose this c corresponds to the relative vorticity. c b(l,m,n) : any scalar field. for our purpose this c corresponds to the streamfunction field. c c dy : grid distance in meters in the c meridional direction. c l : number of grid points in the zonal c direction. c m : number of grid points in the meridional c direction. c n : number of vertical levels. c l1 : l-1 c m1 : m-1 c l2 : l-2 c m2 : m-2 c ladd : number of grid points added to the c eastern boundary to make the domain c cyclic in the zonal direction. c real a(l,m,n), b(l,m,n), dx(m) k = 1 do 10150 j= 2, m1 jp1 = j+1 jm1 = j-1 do 10150 i= 2, l1 10150 a(i,j,k) = ((b(i+1,j,k)+b(i-1,j,k)-2.*b(i,j,k))/dx(j)**2 & +(b(i,jp1,k)+b(i,jm1,k)-2.*b(i,j,k))/dy**2) do 10151 j= 2, m1 a(1,j,k) = ((b(2,j,k)+b(l1,j,k)-2.*b(1,j,k))/dx(j)**2 & +(b(1,j+1,k)+b(1,j-1,k)-2.*b(1,j,k))/dy**2) 10151 a(l,j,k) = a(1,j,k) c c Subroutine BOUND is used to set the southern c and northern boundary conditions. c call BOUND (a,l,m,n,l1,m1,m2,ladd) c c The next do loop is for non-cyclic case only. c if (ladd.eq.0) then do 10152 j = 1, m a(1,j,k) = a(2,j,k) 10152 a(l ,j,k) = a(l1,j,k) endif return end subroutine BOUND (a,l,m,n,l1,m1,m2,ladd) c c This subroutine is used to obtain the southern c and northern boundary values by extrapolation. for c domain that is cyclic in the zonal direction, the c values are extrapolated linearly outwards . for c domain that have zonal boundary values formed by c extending the values 1 grid point outwards, a c similar extension is done here for the southern and c northern boundary values. c c Definitions : c a(l,m,n) : variable whose southern and northern c boundary values have to be obtained by c extrapolation. c l : number of grid points in the zonal c direction. c m : number of grid points in the meridional c direction. c n : number of levels in the vertical. c l1 : l-1 c m1 : m-1 c m2 : m-2 c ladd : number of grid points added to the c eastern boundary to make the domain c cyclic in the zonal direction. c real a(l,m,n) k = 1 c do 10110 i= 1, l a(i,1,k) = 2.*a(i,2,k)-a(i,3,k) a(i,m,k) = 2.*a(i,m1,k)-a(i,m2,k) 10110 continue c c The following loop is for the non-cyclic case only. c if (ladd.eq.0) then do 10111 i= 1, l a(i,1,k) = a(i,2,k) 10111 a(i,m,k) = a(l1,m1,k) endif return end