program box_model c c JIM POTEMRA , University of Hawaii at Manoa c c This program is a 1-1/2 layer model originally developed c for the URI workshop. c No entrainment/detrainment c No thermodynamics c C-Grid used, with j = 0 along bottom row c j = ( jlat + 1 ) / 2 at the equator c j = jlat + 1 along top row c only v's along top and bottom row c i = 0 along first (left) column c i = ilon + 1 along last (right) column c only u's along first and last columns c u(i,j) to the right of h(i,j) c v(i,j) below h(i,j) c For the GrADS display, only the (1,ilon:1,jlat) are saved. c ------------------------- Declare Variables ------------------------- c iflag_force = 0 : No wind stress applied c = 1 : Wind stress applied along western side c iflag_bump = 0 : Flat height field initially c = 1 : Gaussian bump in height field initially c set the length in days of model run c set number of longitude and latitude points [note, jlat should be c an odd number, and the equator will be at (jlat+1)/2] c set time step in minutes c set initial upper layer thickness in meters c set speed of first baroclinic mode Kelvin wave in meters per second c set viscosity in per meters per second c set delta-x in meters c set delta-y in meters c set the number of labels for output color scale c set the min and max for output c set the output resolution (lon [deg], lat [deg], time [day]) c c set center latitude c c --------------------------------------------------------------------- c USER DEFINED PARAMETERS c --------------------------------------------------------------------- parameter ( iflag_force = 1 ) parameter ( iflag_bump = 0 ) parameter ( maxday = 10 ) parameter ( ilon = 319, jlat = 57 ) parameter ( dt_min = 35.0 ) parameter ( h_init = 75.0 ) parameter ( cspd = 2.4 ) parameter ( visc = 0.0 ) parameter ( ddx = 27500.0 ) parameter ( ddy = 27500.0 ) parameter ( labels = 7 ) parameter ( ix_disp = 2, iy_disp = 2, out_int = 1.0 ) parameter ( y_0 = -10.e5 ) parameter ( period_day = 300. ) c parameter ( y_0 = 0. ) c --------------------------------------------------------------------- integer i, j, k, nw, nc, nd, nn, nstep, nt, icount, jcount, + int_label, irecl, irec, orec integer scale(labels) integer mask_u(0:ilon+1,0:jlat), mask_v(1:ilon+1,0:jlat+1) real pi, g, Rearth, omega, dt, beta, y, out_inc, max_damp, + y_u, y_v, u_p1, v_p1, dx_h, dy_h, x_d1, y_d1, eta real u(0:ilon+1,0:jlat,3), v(1:ilon+1,0:jlat+1,3), + u_c(0:ilon+1,0:jlat), v_c(1:ilon+1,0:jlat+1), + u_p(0:ilon+1,0:jlat), v_p(1:ilon+1,0:jlat+1), + u_d(0:ilon+1,0:jlat), v_d(1:ilon+1,0:jlat+1), + u_c1(0:jlat+1), v_c1(0:jlat+1), h(1:ilon+1,0:jlat,3), + wx(0:ilon+1,0:jlat), wy(1:ilon+1,0:jlat+1), & damp(0:jlat+1),damp_x(0:ilon+1) c -------------------------- Open Data Sets --------------------------- irecl = ilon * jlat c open ( 16, file = "bump.dat", form = "unformatted", status='unknown' ) c --------------------------- Set Constants --------------------------- c pi = 4.0 * atan ( 1.0 ) /* pi */ c g = cspd * cspd / h_init /* reduced gravity */ c Rearth = 6378400.0 /* radius of the earth */ c omega = 2.0 * pi / ( 24.0 * 3600.0 ) /* rotation rate of the earth */ c dt = dt_min * 60.0 /* time step in seconds */ c out_inc = out_int * 24.0 * 3600.0 / dt /* output increment */ c beta = 2.0 * omega / Rearth /* beta (for beta-plane) */ c max_damp= 1.0 / ( 24.0 * 3600.0 ) /* max damp for open n/s bound*/ pi = 4.0 * atan ( 1.0 ) g = cspd * cspd / h_init Rearth = 6378400.0 omega = 2.0 * pi / ( 24.0 * 3600.0 ) dt = dt_min * 60.0 out_inc = out_int * 24.0 * 3600.0 / dt beta = 2.0 * omega / Rearth max_damp= 1.0 / ( 24.0 * 3600.0 ) period = period_day*24.*3600. c ----------------------- Initialize All Fields ----------------------- irec = 1 orec = 0 do 10 k = 1, 3 do 20 j = 0, jlat do 20 i = 1, ilon + 1 u(i,j,k) = 0.0 v(i,j,k) = 0.0 h(i,j,k) = h_init 20 continue do 30 i = 1, ilon + 1 v(i,jlat+1,k) = 0.0 30 continue do 40 j = 0, jlat u(0,j,k) = 0.0 40 continue 10 continue do j = 0, jlat do i = 1, ilon + 1 mask_u(i,j) = 1 mask_v(i,j) = 1 enddo mask_u(0,j) = 1 enddo do i = 1, ilon + 1 mask_v(i,jlat+1) = 1 enddo c ------------------- Prescribing the Forcing Terms ------------------- do 50 j = 0, jlat do 50 i = 1, ilon + 1 wx(i,j) = 0.0 wy(i,j) = 0.0 50 continue do 60 i = 1, ilon + 1 wy(i,jlat+1) = 0.0 60 continue do 70 j = 0, jlat wx(0,j) = 0.0 70 continue if ( iflag_force.eq.1 ) then do 80 j = 0, jlat/2 ff = 1. c ff = sin(pi*real(j)/real(jlat)) do 80 i = 162, 238 wx(i,j) = ff* ( -2.5E-5 ) + * (sin ( pi * ( real ( i ) - 162.0 ) / 75. )) wy(i,j) = 0.0 80 continue endif do j = jlat/2+1, jlat+1 do i = 161, 250 mask_v(i,j) =0 enddo enddo do j = jlat/2+1, jlat do i = 160, 250 mask_u(i,j) =0 enddo enddo write(18, '(5F10.4)') + ( ( 1.0*mask_u(i,j), i = 1, ilon ), j = 1, jlat ) write(18, '(5F10.4)') + ( ( 1.e5*wx(i,j), i = 1, ilon ), j = 1, jlat ) c ----------------------- Set Intial Conditions ----------------------- if ( iflag_bump.eq.1 ) then do 90 j = 1, jlat do 90 i = 100, 160 y = 2.0 * ddy + * ( real ( j ) - real ( jlat + 1 ) / 2.0 ) h(i,j,1) = 20.0 * exp ( -beta * y * y / ( 2.0 * cspd ) ) + * sin ( ( real ( i ) - 100.0 ) * pi / 60.0 ) + + 200.0 h(i,j,2) = h(i,j,1) 90 continue endif c ----------- Pre-compute Some Frequently Used Parameters ------------- do 100 j = 0, jlat + 1 y_u = 2.0 * ddy + * ( real ( j ) - real ( jlat ) / 2.0 ) + y_0 y_v = 2.0 * ddy + * ( real ( j ) - real ( jlat + 1 ) / 2.0 ) + y_0 u_c1(j) = beta * y_u v_c1(j) = -beta * y_v 100 continue u_p1 = -g / ( 2.0 * ddx ) v_p1 = -g / ( 2.0 * ddy ) dx_h = 1.0 / ( 2.0 * ddx ) dy_h = 1.0 / ( 2.0 * ddy ) x_d1 = visc / ( 4.0 * ddx * ddx ) y_d1 = visc / ( 4.0 * ddy * ddy ) c ------------- Set Damping for North South Open Boundaries ----------- do 110 j = 5, jlat - 5 damp(j) = 0.0 110 continue do 120 j = 0, 4 damp(j) = -max_damp * real ( 5 - j ) / 5.0 damp(j+jlat-3)= -max_damp * real ( j+1 ) / 5.0 120 continue damp(jlat+1)= -max_damp do i = 5, ilon - 5 damp_x(i) = 0.0 enddo do i = 0, 4 damp_x(i) = -max_damp * real ( 5 - i ) / 5.0 damp_x(i+ilon-3)= -max_damp * real ( i+1 ) / 5.0 enddo c ------------------ Initialize for Time Integration ------------------ nw = 1 nc = 2 nd = 3 nn = 0 nstep = 0 tme = 0. c ---------------------- Begin Time Integration ----------------------- 130 nn = nn + 1 nstep = 0 140 nstep = nstep + 1 nt = nw nw = nd nd = nc nc = nt tme = tme + dt c ---------------- compute wind forcing factor ----------------------- w_fac = sin(2.*pi*tme/period) c ---------------- Compute the Pressure Gradient Terms ---------------- do 150 j = 0, jlat u_p(0,j) = 0.0 u_p(ilon+1,j) = 0.0 do 150 i = 1, ilon u_p(i,j) = u_p1 * ( h(i+1,j,nc) - h(i,j,nc) ) * h_init 150 continue do 160 i = 1, ilon + 1 v_p(i,0) = 0.0 v_p(i,jlat+1) = 0.0 do 160 j = 1, jlat v_p(i,j) = v_p1 * ( h(i,j,nc) - h(i,j-1,nc) ) * h_init 160 continue c -------------------- Compute the Coriolis Terms --------------------- do 170 j = 0, jlat u_c(0,j) = u_c1(j) * 0.5 * ( v(1,j,nc) + v(1,j+1,nc) ) u_c(ilon+1,j) = u_c1(j) * 0.5 + * ( v(ilon+1,j,nc) + v(ilon+1,j+1,nc) ) do 170 i = 1, ilon u_c(i,j) = u_c1(j) * 0.25 + * ( v(i,j,nc) + v(i+1,j,nc) + + v(i+1,j+1,nc) + v(i,j+1,nc) ) 170 continue u_c(0,0) = 0.0 do 180 i = 1, ilon + 1 c v_c(i,0) = v_c1(0) * 0.5 * ( u(i-1,0,nc) + u(i,0,nc) ) c v_c(i,jlat+1) = v_c1(jlat+1) * 0.5 c + * ( u(i-1,jlat,nc) + u(i,jlat,nc) ) v_c(i,0) = - 0.5 * u_c1(0)* ( u(i-1,0,nc) + u(i,0,nc) ) v_c(i,jlat+1) = - u_c1(jlat) * 0.5 + * ( u(i-1,jlat,nc) + u(i,jlat,nc) ) do 180 j = 1, jlat c v_c(i,j) = v_c1(j) * 0.25 c + * ( u(i-1,j-1,nc) + u(i,j-1,nc) c + + u(i,j,nc) + u(i-1,j,nc) ) v_c(i,j) = - 0.25 + * ( (u(i-1,j-1,nc) + u(i,j-1,nc))*u_c1(j-1) + + (u(i,j,nc) + u(i-1,j,nc))*u_c1(j) ) 180 continue v_c(ilon+1,jlat+1) = 0.0 c -------------------- Compute the Viscous Terms ---------------------- c c do 190 j = 1, jlat - 1 c do 190 i = 1, ilon c u_d(i,j) = x_d1 c + * ( u(i+1,j,nd) + u(i-1,j,nd) - 2.0 * u(i,j,nd) ) c + + y_d1 c + * ( u(i,j+1,nd) + u(i,j-1,nd) - 2.0 * u(i,j,nd) ) c190 continue c do 200 j = 1, jlat c do 200 i = 2, ilon c v_d(i,j) = x_d1 c + * ( v1(i+1,j,nd) + v1(i-1,j,nd) - 2.0 * v1(i,j,nd) ) c + + y_d1 c + * ( v1(i,j+1,nd) + v1(i,j-1,nd) - 2.0 * v1(i,j,nd) ) c200 continue c do 210 j = 0, jlat c u_d(0,j) = c u_d(ilon+1,j) = c v_d(1,j) = c v_d(ilon+1,j) = c210 continue c do 220 i = 1, ilon + 1 c u_d(i,jlat+1) = c u_d(i,0) = c v_d(i,jlat+1) = c v_d(i,0) = c220 continue c u_d(0,jlat+1) = c u_d(0,0) = c v_d(ilon+1,jlat+1) = c v_d(1,jlat+1) = c c ---------------------- Update the u,v,h fields ---------------------- do 230 j = 0, jlat do 230 i = 0, ilon + 1 u(i,j,nw) = u(i,j,nd) + 2.0 * dt * mask_u(i,j) + * ( u_c(i,j) + u_p(i,j) + u_d(i,j) + w_fac*wx(i,j) + + (damp(j)+damp_x(i)) * u(i,j,nd) ) 230 continue do 240 j = 0, jlat + 1 do 240 i = 1, ilon + 1 v(i,j,nw) = v(i,j,nd) + 2.0 * dt * mask_v(i,j) + * ( v_c(i,j) + v_p(i,j) + v_d(i,j) + w_fac*wy(i,j) + + (damp(j)+damp_x(i)) * v(i,j,nd) ) 240 continue do 250 j = 0, jlat do 250 i = 1, ilon + 1 h(i,j,nw) = h(i,j,nd) + - 2.0 * dt * dx_h * ( u(i,j,nc) - u(i-1,j,nc) ) + - 2.0 * dt * dy_h * ( v(i,j+1,nc) - v(i,j,nc) ) + + 2.0 * dt * (damp(j)+damp_x(i))*(h(i,j,nd)-h_init) 250 continue c ------------- Do Next Time Step, Unless Time for Output ------------- if ( nstep.lt.int(out_inc) ) go to 140 c ------------------------- Set Output Values ------------------------- do 260 j = 0, jlat do 260 i = 0, ilon + 1 u(i,j,nc) = 0.5 * ( u(i,j,nw) + u(i,j,nc) ) u(i,j,nw) = u(i,j,nc) 260 continue do 270 j = 0, jlat + 1 do 270 i = 1, ilon + 1 v(i,j,nw) = 0.5 * ( v(i,j,nw) + v(i,j,nc) ) v(i,j,nc) = v(i,j,nw) 270 continue do 280 j = 0, jlat do 280 i = 1, ilon + 1 h(i,j,nc) = 0.5 * ( h(i,j,nw) + h(i,j,nc) ) h(i,j,nw) = h(i,j,nc) 280 continue print *,'Computed Model Day ', nn c --------------------- Output Results (GRADS Format) ------------------ if ( mod(nn,1).eq.0.or.nn.eq.1 ) then orec = orec + 1 write ( 16, '(5F10.4)' ) + ( ( u(i,j,nw), i = 1, ilon ), j = 1, jlat ) write ( 17, '(5F10.4)' ) + ( ( v(i,j,nw), i = 1, ilon ), j = 1, jlat ) write ( 18, '(5F10.4)' ) + ( ( h(i,j,nw), i = 1, ilon ), j = 1, jlat ) write(6,*) orec c write(6,*) ( ( h(i,j,nw), i = 1, ilon ), j = 1, jlat ) endif c ---------------- Do Next Time Step, Unless End of Run ---------------- if ( nn.lt.maxday ) go to 130 print *,"Number of times written:", orec stop end