c This program is an atmospheric or meteorological simulator that c simulates the transport of momentum, heat and potentially other c quantities in a stably to neutrally stratified atmosphere. It uses c basically conservative advection/diffusion method for tracking the c different physical quantities. c The Task list for future development includes: c (1) Separation of vertical and horizontal diffusion with c fine time steps for z-diffusion c (2) Include velocity criterion for neutrally stratified c mixing. c (3) Do output for calibration stations so that parameters c can be calibrated c (4) Automate parameter calibration Program heatcalc character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 c write(*,*) "Can I get started?" c Initialize the different model values including parameter c data inputs, and the initial temperature, heat, and velocity c fields. call initvalues c--------Loop Through Different Cases c For the most exterior loop we go through the different c parameter cases in search of the ideal paramters for the c simulation model ncmx=4 ndaymx=5 do 5 nc=2,ncmx call initcase(nc) c--------Simulation Loop with Printing: c Now we do the simulation loop, where we go through stepping c time and periodically outputting data do 10 it=1,itmx c Here we print out the results periodically to an output file c these can then be picked up by another program and used for c display and analysis. if (time.gt.72) then call printvalues endif c -------Core Simulation Loop: do 15 nd=1,ndit c ....Calculation of vertical velocities call vzcalc c ....Calculation of Pressure gradients etc. call momdrivers c ....Heat advection subroutine call heatadvect c ....Advection of momentum density call thadvect c ....Advection of momentum density call momadvect c ....Acceleration of the air flow call accelerate c ....Heat sources and inputs call heatinputs c ....Momentum Diffusion Routine call momdiffuse write(*,*) "time =",time c ....Bulk Radiative Heat loss Routine call radloss c ....Heat Diffusion Routine call heatdiffuse time=time+dt/3600 c jump out of the loop if we hit the maximum time. if (time/24.gt.ndaymx) goto 5 write(*,*) "time =",time 15 continue c ------- End Core Simulation Loop: 10 continue c-------- End Simulation Loop with Printing: 5 continue c-------- End Parameter cases loop stop end c ------------------End Main Program --------------------- c------------------Start function realz ----------- function realz(z) zmax=9000 realz = -zmax*log(1-(z/zmax)) return end c------------------End function realz ------------- c------------------Start function densz ----------- function densz(z) zmax=9000 densz = zmax*(1-exp(-z/zmax)) return end c------------------End function realz ------------- c------------------Start subroutine initvalues ----------- subroutine initvalues character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 character*1 iz(1440000) real izvalav c Units for parameters c c Please note that the units and descriptions of the different c parameters are as follows: c c time is the current time in hours in local time c dx, dy, and dz are the obvious spacial grid spacings in c meters, except that note the z-coordinate is in what I call c density height. It is a z-interval that is inversely c proportional to the mean air density at that height. c Tup, Tsea, Tz0 are temperature values in degrees C where c Tup is the mean value of the temperature at the upper c boundary, dTup is the magnitude of the daily variation c of the upper boundary temperature, and Tz0 is the initial c temperature of the ground. c ittotal, itmax, ndit ittotal is the total number of iterations c of the simulation when it is completed (a bit less than 10 days c is probably enough, while itmx is the number of printing c iterations, and ndit is the number of non-printing iterations c between printing iterations. c write(*,*) " in initval" zmax=9000. time=6. ittot=100000 ndit=16 itmx=ittot/ndit nxsea=40 nxcool=7 nxmx=60 nymx=60 nzmx=40 dz=6000/nzmx nx0=1 ny0=1 nz0=1 pi=4.*atan(1.) dTsea=14. Tsea=27.-dTsea*0.5 dTatm=6. Vxup=0. Vyup=0. Tz0=14. Tcool=22. Z0=-600. ratelaps= -0.01 dx=1200000/nxmx dy=1200000/nymx Zup=Z0+dz*(nzmx-1) st0=-0.0060 Tup=Tz0+dTatm+st0*realz(Zup) dt=900. area=dx*dy c Now we reset dt to reflect the requirements of the advection c step. dspace=sqrt(area) vmax=40. ndt=1+int((vmax*dt)/dspace) dt=dt/ndt c Density in kilograms per cubic meter at sea level Dens0 = 1.225 c values of the order of a few ten thousandths c of a pascal per meter are about right. GPx0=0.0000 GPy0=0.0000 c Next convert to gradients in stratification c dt0x=-GPx0*300/(Dens0*9.8*dz*nzmx) c dt0y=-GPy0*300/(Dens0*9.8*dz*nzmx) dt0x=0. dt0y=0. Czmx= 15. Czmn= 2. c For Czgd, we are going to use it as a conductivity c coefficient for the near-ground layer. We multiply c this coefficient, times the heat transport rate out c of the ground and multiply by the thickness of the c ground layer (zgd) to get the temperature difference between c the ground layer and the neutral layer. c The heat transport rate out of the ground is approximated c by Czgd=0.1 dzgd=10. Cxy = 5000. vmix=5. Cvzmx= 10. Cvzmn= 1.5 c Cvzgd= 0.008 Cvzgd= 0.02 Cvzsea= 0.0005 Cvxy1 = 10000. Cvxy2 = 3000. c Heat capacity in Joules per kilogram per degree Cp = 1004. Cpsea = 10000. Rad0 = 1000. c Rnt is the heat loss rate due to radiation to space c assume that it is about 10% of the incoming solar radiation c Rnt = -25. c This version of Rnt is the factor that when multiplied times c the difference between the Temperature and the equilibrium c temperature gives the outgoing heat radiation Rnt=Cp*Dens0*(200)/(4*3600) abland = 0.7 absea = 0.2 Cpfact=Dens0*Cp*area*dz Cpfacti=1/Cpfact alpha=0.2 do 10 i=1,nxmx do 20 j=1,nymx c at deg degrees north latitude... the f factor deg=15.+10.*float(j-nymx/2)/float(nymx) fc(j)=2.*sin(deg*pi/180)/(24*3600) c dt0(i,j)=float(i-nxmx/2)*dx*dt0x+float(j-nymx/2)*dy*dt0y dt0(i,j)=float(i-nxmx/2)*dt0x/nxmx+float(j-nymx/2)*dt0y/nymx do 30 k=1,nzmx elev=float(k)*dz+Z0 T(i,j,k)=Tz0+dt0(i,j)+dTatm+st0*realz(elev) Q(i,j,k)=Cpfact*T(i,j,k) vx(i,j,k)=0. vy(i,j,k)=0. vz(i,j,k)=0. grdpx(i,j,k)=0. grdpy(i,j,k)=0. rlp(k)=ratelaps/(1-elev/zmax) 30 continue vx(i,j,nzmx)=Vxup vy(i,j,nzmx)=Vyup th(i,j)=0. 20 continue 10 continue c Now we try to read in the digital elevation data. ndata=1440000 113 format(1440000a1) open(unit=13,file="e3545n1020.dem") read(13,113) (iz(m),m=1,ndata) close(unit=13) ksea=-Z0/dz+1 ndx=1200/nxmx ndy=1200/nymx ndtot=ndx*ndy do 40 i=1,nxmx do 50 j=1,nymx izvalav=0 nizval=0 do 41 nx=1,ndx do 51 ny=1,ndy m=((1200*i)/nxmx+nx-1)+1200*(1200-ny+1-((1200*j)/nymx)) c izval=int(iz(m)) izval=ichar(iz(m)) if (izval.le.250) then izvalav=izvalav+izval nizval=nizval+1 endif 51 continue 41 continue if (nizval.gt.(ndtot/2)) then izvalav=izvalav/(nizval) else izvalav=254 endif c m=((1200*i)/nxmx+ndx/2-1)+1200*(1200-ndy+1-((1200*j)/nymx)) c izval=int(iz(m)) if (izvalav.gt.250) then kmn(i,j)=ksea dep(i,j)=0. nsea(i,j)=1 else dep(i,j)=densz(20*(izvalav-30))-Z0 kmn(i,j)=(dep(i,j)/dz) + 1 dep(i,j)=kmn(i,j)*dz-dep(i,j) nsea(i,j)=0 endif ntmp(i,j)=0 kmin=kmn(i,j) elev=dz*kmin+Z0 T(i,j,kmin)=T(i,j,kmin)-st0*dep(i,j)/(1-elev/zmax) Q(i,j,kmin)=Cpfact*T(i,j,kmin) 50 continue 40 continue c Now lets do side boundaries do 45 j=1,nymx kmn(1,j)=kmn(2,j) nsea(1,j)=nsea(2,j) ntmp(1,j)=ntmp(2,j) dep(1,j)=dep(2,j) kmn(nxmx,j)=kmn(nxmx-1,j) nsea(nxmx,j)=nsea(nxmx-1,j) ntmp(nxmx,j)=ntmp(nxmx-1,j) dep(nxmx,j)=dep(nxmx-1,j) 45 continue do 55 i=1,nxmx kmn(i,1)=kmn(i,2) nsea(i,1)=nsea(i,2) ntmp(i,1)=ntmp(i,2) dep(i,1)=dep(i,1) kmn(i,nymx)=kmn(i,nymx-1) nsea(i,nymx)=nsea(i,nymx-1) ntmp(i,nymx)=ntmp(i,nymx-1) dep(i,nymx)=dep(i,nymx-1) 55 continue c calculate topographic slope values dxi2=0.5/dx dyi2=0.5/dy do 46 i=2,nxmx-1 do 56 j=2,nymx-1 slpx(i,j)=(dz*(kmn(i+1,j)-kmn(i-1,j)) 1 -dep(i+1,j)+dep(i-1,j))*dxi2 slpy(i,j)=(dz*(kmn(i,j+1)-kmn(i,j-1)) 1 -dep(i,j+1)+dep(i,j-1))*dyi2 c Do geostrophic velocities kminp=kmn(i,j)+1 do 66 k=kminp,nzmx vx(i,j,k)=-GPy0/(fc(j)*Dens0) vy(i,j,k)=GPx0/(fc(j)*Dens0) 66 continue 56 continue 46 continue do 47 j=1,nymx slpx(1,j)=0. slpy(1,j)=slpy(2,j) slpx(nxmx,j)=0. slpy(nxmx,j)=slpy(nxmx-1,j) 47 continue do 57 i=1,nxmx slpx(i,1)=slpx(i,2) slpy(i,1)=0. slpx(i,nymx)=slpx(i,nymx-1) slpy(i,nymx)=0. 57 continue c Do corner slopes: slpy(nxmx,nymx)=0. slpy(1,nymx)=0. slpy(nxmx,1)=0. slpy(1,1)=0. slpx(nxmx,nymx)=0. slpx(1,nymx)=0. slpx(nxmx,1)=0. slpx(1,1)=0. return end c------------------End subroutine initvalues ----------- c------------------Start subroutine printvalues ----------- c c This subroutine is supposed to output some of the data in c a form that can be used by other programs for further c processing. subroutine printvalues character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 common/refined/ref(490,490),tref(490,490) character*1 d1,d2,d3 character*3 ctime c Define arrays for averages dimension vav(490,490,3),vxav(490,490,3),vyav(490,490,3) c Derive the character expression for the time nc1=48+int(time/100) nc2=48+int(time/10)-10*int(time/100) nc3=48+int(time)-10*(nc2-48)-100*(nc1-48) d1=char(nc1) d2=char(nc2) d3=char(nc3) ctime=d1//d2//d3 c We have a running average factor which we set here avfact=0.03 c Output the different quantitites for the different reference c locations 111 format(a10,a5,f8.3,a6,f8.3,a6,f8.3,a6,f8.0,a5,f8.3) c Output for Asmara i=(38.9167-35)*120000/dx j=(15.28333-10)*120000/dy kmin=kmn(i,j) elev=realz((kmin)*dz+Z0) elev=elev-dep(i,j)/(1-elev/zmax) write(11,111) "Asmera: "," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time write(*,111) "Asmera: "," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time 119 format(a5,f6.2,a6,f6.2,a6,f6.2,a6,f6.2,a6,f6.0,a5,f8.3) c Output verticle profile above the station open(unit=10,file="Asm"//ctime//".txt") do 901 k=kmin,nzmx elev=realz((k)*dz+Z0) pTemp=T(i,j,k)+0.01*elev write(10,119) " T = ",T(i,j,k) 1 ," pT = ",pTemp 2 ," Vx = ",vx(i,j,k), " Vy = ",vy(i,j,k)," El = ",elev 3 ," t = ",time 901 continue close(unit=10) c Output for Massawa i=(39.43-35)*120000/dx j=(15.6-10)*120000/dy kmin=kmn(i,j) elev=realz((kmin)*dz+Z0) elev=elev-dep(i,j)/(1-elev/zmax) write(11,111) "Massawa: "," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time write(*,111) "Massawa: "," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time c Output verticle profile above the station open(unit=10,file="Mas"//ctime//".txt") elev=realz((kmin)*dz+Z0) pTemp=Tg(i,j)+0.01*elev c write(10,119) " T = ",Tg(i,j) c 1 ," pT = ",pTemp c 2 ," Vx = ",vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev c 3 ," t = ",time do 904 k=kmin,nzmx elev=realz((k)*dz+Z0) pTemp=T(i,j,k)+0.01*elev write(10,119) " T = ",T(i,j,k) 1 ," pT = ",pTemp 2 ," Vx = ",vx(i,j,k), " Vy = ",vy(i,j,k)," El = ",elev 3 ," t = ",time 904 continue close(unit=10) c Output for Aseb i=(42.7-35)*120000/dx j=(13.1-10)*120000/dy kmin=kmn(i,j) elev=realz((kmin)*dz+Z0) elev=elev-dep(i,j)/(1-elev/zmax) write(11,111) "Aseb: "," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time write(*,111) "Aseb: "," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time c Output verticle profile above the station open(unit=10,file="Asb"//ctime//".txt") elev=realz((kmin)*dz+Z0) pTemp=Tg(i,j)+0.01*elev c write(10,119) " T = ",Tg(i,j) c 1 ," pT = ",pTemp c 2 ," Vx = ",vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev c 3 ," t = ",time do 903 k=kmin,nzmx elev=realz((k)*dz+Z0) pTemp=T(i,j,k)+0.01*elev write(10,119) " T = ",T(i,j,k) 1 ," pT = ",pTemp 2 ," Vx = ",vx(i,j,k), " Vy = ",vy(i,j,k)," El = ",elev 3 ," t = ",time 903 continue close(unit=10) c Output for South Red Sea i=(43.0-35)*120000/dx j=(13.3-10)*120000/dy kmin=kmn(i,j) elev=realz((kmin)*dz+Z0) elev=elev-dep(i,j)/(1-elev/zmax) write(11,111) "S Red Sea:"," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time write(*,111) "S Red Sea:"," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time c Output verticle profile above the station open(unit=10,file="SRS"//ctime//".txt") elev=realz((kmin)*dz+Z0) pTemp=Tg(i,j)+0.01*elev do 905 k=kmin,nzmx elev=realz((k)*dz+Z0) pTemp=T(i,j,k)+0.01*elev write(10,119) " T = ",T(i,j,k) 1 ," pT = ",pTemp 2 ," Vx = ",vx(i,j,k), " Vy = ",vy(i,j,k)," El = ",elev 3 ," t = ",time 905 continue close(unit=10) c Output for Tessenei i=(36.7-35)*120000/dx j=(15.06-10)*120000/dy kmin=kmn(i,j) elev=realz((kmin)*dz+Z0) elev=elev-dep(i,j)/(1-elev/zmax) write(11,111) "Tessenei: "," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time write(*,111) "Tessenei: "," T = ",T(i,j,kmin)," Vx = ", 1 vx(i,j,kmin), " Vy = ",vy(i,j,kmin)," El = ",elev 2 ," t = ",time c Output verticle profile above the station open(unit=10,file="Tes"//ctime//".txt") elev=realz((kmin)*dz+Z0) pTemp=Tg(i,j)+0.01*elev c write(10,119) " T = ",Tg(i,j) c 1 ," pT = ",pTemp c 2 ," Vx = ",vx(i,j,k), " Vy = ",vy(i,j,k)," El = ",elev c 3 ," t = ",time do 902 k=kmin,nzmx elev=realz((k)*dz+Z0) pTemp=T(i,j,k)+0.01*elev write(10,119) " T = ",T(i,j,k) 1 ," pT = ",pTemp 2 ," Vx = ",vx(i,j,k), " Vy = ",vy(i,j,k)," El = ",elev 3 ," t = ",time 902 continue close(unit=10) c Let's output the temperature, and potential temperature c We will do this now as one=byte pixel images. 222 format(a1,$) 223 format(a2) 224 format(i3) c Output the time and other relevant parameters open(unit=10,file="temp"//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx,nzmx write(10,224) 255 do 10 k=nzmx,1,-1 c do 20 j=nymx,1,-1 j=nymx/2 do 21 i=1,nxmx write(10,222) int(4*T(i,j,k)+100) 21 continue 20 continue 10 continue close(unit=10) c Output the time and other relevant parameters open(unit=10,file="ptemp"//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx,nzmx write(10,224) 255 do 40 k=nzmx,1,-1 c do 50 j=nymx,1,-1 j=nymx/2 do 60 i=1,nxmx tmp(i,j,k)=T(i,j,k)+0.01*realz(((k-1)*dz+Z0)) 60 continue do 51 i=1,nxmx write(10,222) int(4*tmp(i,j,k)+0) 51 continue 50 continue 40 continue close(unit=10) c Output the x-component of the velocity open(unit=10,file="velx"//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx,nzmx write(10,224) 255 do 70 k=nzmx,1,-1 c do 80 j=nymx,1,-1 j=nymx/2 do 81 i=1,nxmx write(10,222) int(5*vx(i,j,k)+100) 81 continue 80 continue 70 continue close(unit=10) c Output the x-component of the pressure gradient open(unit=10,file="gradpx"//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx,nzmx write(10,224) 255 do 90 k=nzmx,1,-1 c do 100 j=nymx,1,-1 j=nymx/2 do 101 i=1,nxmx write(10,222) int(10000*grdpx(i,j,k)+100) 101 continue 100 continue 90 continue close(unit=10) c Output the z-component of the velocity open(unit=10,file="velz"//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx,nzmx write(10,224) 255 do 110 k=nzmx,1,-1 c do 120 j=nymx,1,-1 j=nymx/2 do 121 i=1,nxmx write(10,222) int(50*vz(i,j,k)+100) 121 continue 120 continue 110 continue close(unit=10) do 99 nup=0,2,2 if (nup.eq.0) then d1="A" else if (nup.eq.2) then d1="B" else d1="C" endif endif c Lets set parameters for refined grid data: nref=3 nxmx2=nxmx*nref nymx2=nymx*nref c We refine the velocity fields before we output them do 131 i=1,nxmx do 141 j=1,nymx k=kmn(i,j)+nup c k=(200+1000*(nup-1)-Z0)/dz tref(i,j)=vx(i,j,k+1)*(1-dep(i,j)/dz)+vx(i,j,k)*dep(i,j)/dz 141 continue 131 continue call refine(ref,tref,nxmx,nymx,490,490,490,490,nref) c Output the x-component of the velocity open(unit=10,file="velxxy"//d1//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 140 j=nymx2,1,-1 do 142 i=1,nxmx2 vxav(i,j,nup+1)=vxav(i,j,nup+1)*(1-avfact)+avfact*ref(i,j) write(10,222) int(5*ref(i,j)+100) 142 continue 140 continue close(unit=10) c Output the average open(unit=10,file="vxaxy"//d1//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 144 j=nymx2,1,-1 do 145 i=1,nxmx2 write(10,222) int(5*vxav(i,j,nup+1)+100) 145 continue 144 continue close(unit=10) c Finished refined VX output c We refine the velocity fields before we output them do 151 i=1,nxmx do 161 j=1,nymx k=kmn(i,j)+nup c k=(200+1000*(nup-1)-Z0)/dz tref(i,j)=vy(i,j,k+1)*(1-dep(i,j)/dz)+vy(i,j,k)*dep(i,j)/dz 161 continue 151 continue nref=3 call refine(ref,tref,nxmx,nymx,490,490,490,490,nref) c Output the y-component of the velocity open(unit=10,file="velyxy"//d1//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 160 j=nymx2,1,-1 do 162 i=1,nxmx2 vyav(i,j,nup+1)=vyav(i,j,nup+1)*(1-avfact)+avfact*ref(i,j) write(10,222) int(5*ref(i,j)+100) 162 continue 160 continue close(unit=10) c Output the magnitude open(unit=10,file="vyaxy"//d1//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 164 j=nymx2,1,-1 do 165 i=1,nxmx2 write(10,222) int(5*vyav(i,j,nup+1)+100) 165 continue 164 continue close(unit=10) c We refine the velocity fields before we output them do 171 i=1,nxmx do 181 j=1,nymx k=kmn(i,j)+nup c k=(200+1000*(nup-1)-Z0)/dz vk=sqrt(vx(i,j,k)**2+vy(i,j,k)**2) vkp=sqrt(vx(i,j,k+1)**2+vy(i,j,k+1)**2) tref(i,j)=vkp*(1-dep(i,j)/dz)+vk*dep(i,j)/dz 181 continue 171 continue nref=3 call refine(ref,tref,nxmx,nymx,490,490,490,490,nref) c Output the magnitude open(unit=10,file="velxy"//d1//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 180 j=nymx2,1,-1 do 182 i=1,nxmx2 vav(i,j,nup+1)=vav(i,j,nup+1)*(1-avfact)+avfact*ref(i,j) write(10,222) int(5*ref(i,j)+50) 182 continue 180 continue close(unit=10) c Output the magnitude open(unit=10,file="vavxy"//d1//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 184 j=nymx2,1,-1 do 185 i=1,nxmx2 write(10,222) int(10*vav(i,j,nup+1)+50) 185 continue 184 continue close(unit=10) c We refine the verticle velocity fields before we output them do 251 i=1,nxmx do 261 j=1,nymx k=nzmx*(nup+2)/5 tref(i,j)=vz(i,j,k) 261 continue 251 continue nref=3 call refine(ref,tref,nxmx,nymx,490,490,490,490,nref) c Output the magnitude open(unit=10,file="vzxy"//d1//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 260 j=nymx2,1,-1 do 262 i=1,nxmx2 write(10,222) int(50*ref(i,j)+100) 262 continue 260 continue close(unit=10) 99 continue c We refine the temperature fields before we output them do 271 i=1,nxmx do 281 j=1,nymx tref(i,j)=T(i,j,kmn(i,j)) 281 continue 271 continue nref=3 call refine(ref,tref,nxmx,nymx,490,490,490,490,nref) c Output the magnitude open(unit=10,file="Txy"//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 280 j=nymx2,1,-1 do 282 i=1,nxmx2 write(10,222) int(2*ref(i,j)+100) 282 continue 280 continue close(unit=10) c Output the magnitude of the potential temperature open(unit=10,file="pTxy"//ctime//".pgm") write(10,223) "P5" write(10,*) nxmx2,nymx2 write(10,224) 255 do 290 j=nymx2,1,-1 do 292 i=1,nxmx2 i2=i/nref+1 j2=j/nref+1 pTemp=ref(i,j)+0.01*realz((dz*(kmn(i2,j2))+Z0-dep(i2,j2))) write(10,222) int(4*pTemp+0) 292 continue 290 continue close(unit=10) return end c------------------End subroutine printvalues ----------- c------------------Start subroutine momdrivers ----------- subroutine momdrivers character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 grav=9.8 cmprs=1./300. pfact=grav*dz*cmprs*Dens0 do 10 i=1,nxmx do 20 j=1,nymx dPres=grav*Dens0*(1+(Tz0-T(i,j,nzmx))*cmprs)*th(i,j) kmin=kmn(i,j) do 30 k=nzmx,kmin,-1 dPres=dPres+(Tz0-T(i,j,k))*pfact tmp(i,j,k)=dPres 30 continue if (kmin.ge.2) then kminm=kmin-1 do 35 k=kminm,1,-1 dPres=dPres+(Tz0-T(i,j,k))*pfact tmp(i,j,k)=dPres 35 continue endif 20 continue 10 continue dxi2=1./(2.*dx) dyi2=1./(2.*dy) do 40 i=2,nxmx-1 do 50 j=2,nymx-1 kmin=kmn(i,j) do 60 k=kmin,nzmx grdpx(i,j,k)=dxi2*(tmp(i+1,j,k)-tmp(i-1,j,k)) grdpy(i,j,k)=dyi2*(tmp(i,j+1,k)-tmp(i,j-1,k)) if ((k.lt.kmn(i+1,j)).or. 1 (k.lt.kmn(i-1,j))) then grdpx(i,j,k)=0. endif if ((k.lt.kmn(i,j+1)).or. 1 (k.lt.kmn(i,j-1))) then grdpy(i,j,k)=0. endif 60 continue 50 continue 40 continue c Now we do the gradients at the boundaries do 70 k=1,nzmx do 80 i=1,nxmx grdpx(i,1,k)=grdpx(i,2,k) grdpx(i,nymx,k)=grdpx(i,nymx-1,k) grdpy(i,1,k)=grdpy(i,2,k) grdpy(i,nymx,k)=grdpy(i,nymx-1,k) 80 continue do 90 j=1,nymx grdpx(1,j,k)=grdpx(2,j,k) grdpx(nxmx,j,k)=grdpx(nxmx-1,j,k) grdpy(1,j,k)=grdpy(2,j,k) grdpy(nxmx,j,k)=grdpy(nxmx-1,j,k) 90 continue 70 continue return end c------------------End subroutine momdrivers ----------- c------------------Start subroutine vzcalc ----------- subroutine vzcalc character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 dzovrdx=dz/dx dzovrdy=dz/dy dtwave=0.5*dt do 10 i=2,nxmx-1 do 20 j=2,nymx-1 vzup=0. kmin=kmn(i,j) c Do slope effects near the ground vzdwn=vx(i,j,kmin)*slpx(i,j)+vy(i,j,kmin)*slpy(i,j) do 30 k=kmin,nzmx vzup=vzdwn-(vx(i+1,j,k)-vx(i-1,j,k))*dzovrdx 1 -(vy(i,j+1,k)-vy(i,j-1,k))*dzovrdy vz(i,j,k)=0.5*(vzup+vzdwn) vzdwn=vzup 30 continue th(i,j)=th(i,j)+vz(i,j,nzmx)*dtwave 20 continue 10 continue c At the boundaries enforce zero derivative of the z-velocity do 40 k=1,nzmx do 50 j=1,nymx vz(1,j,k)=vz(2,j,k) vz(nxmx,j,k)=vz(nxmx-1,j,k) if (k.eq.1) then th(1,j)=th(2,j) th(nxmx,j)=th(nxmx-1,j) endif 50 continue do 60 i=1,nxmx vz(i,1,k)=vz(i,2,k) vz(i,nymx,k)=vz(i,nymx-1,k) if (k.eq.1) then th(i,1)=th(i,2) th(i,nymx)=th(i,nymx-1) endif 60 continue 40 continue c Now the following loops are for diffusing the top surface c height displacement. Because really we want it to solve c a Laplace's equation, so if the curvature equals the upwelling, c this should be about right. nsmooth=10 do 99 n=1,nsmooth do 70 i=2,nxmx-1 do 80 j=2,nymx-1 tmp(i,j,1)=0.2*(th(i,j)+th(i+1,j)+th(i-1,j)+th(i,j+1) 1 +th(i,j-1)) 80 continue 70 continue do 71 i=2,nxmx-1 j=1 c tmp(i,j,1)=tmp(i,j+1,1) tmp(i,j,1)=0. j=nymx c tmp(i,j,1)=tmp(i,j-1,1) tmp(i,j,1)=0. 71 continue do 81 j=1,nymx i=1 c tmp(i,j,1)=tmp(i+1,j,1) tmp(i,j,1)=0. i=nxmx c tmp(i,j,1)=tmp(i-1,j,1) tmp(i,j,1)=0. 81 continue do 72 i=1,nxmx do 82 j=1,nymx th(i,j)=tmp(i,j,1) 82 continue 72 continue 99 continue return end c------------------End subroutine vzcalc ----------- c------------------Start subroutine heatadvect ----------- subroutine heatadvect character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 case=1. call advectxy(Q,vx,vy,tmp,nxmx,nymx,nzmx,245,245,95,dx,dy,dt, 1 case) do 40 i=1,nxmx do 50 j=1,nymx kmin=kmn(i,j) do 60 k=kmin,nzmx T(i,j,k)=Cpfacti*Q(i,j,k) 60 continue 50 continue 40 continue call Qadvectz(Q,T,vz,tmp,nxmx,nymx,nzmx,245,245,95,dz,dt) return end c------------------End subroutine heatadvect ----------- c------------------Begin Qadvectz subroutine -------------- subroutine Qadvectz(Q,T,vz,tmp,nxmx,nymx,nzmx,nxd,nyd,nzd, 1 dz,dt) dimension Q(nxd,nyd,nzd),T(nxd,nyd,nzd),vz(nxd,nyd,nzd) dimension tmp(nxd,nyd,nzd) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) dzi=1/dz do 10 i=1,nxmx do 20 j=1,nymx c T(i,j,nzmx)=1.25*T(i,j,nzmx-1)-0.25*T(i,j,nzmx-5) T(i,j,nzmx)=Tup Q(i,j,nzmx)=Cpfact*T(i,j,nzmx) kminp=kmn(i,j)+1 do 30 k=kminp,nzmx-1 dpTdzp=-0.5*(rlp(k+1)+rlp(k))*dz dpTdzm=-0.5*(rlp(k)+rlp(k-1))*dz vzdt=vz(i,j,k)*dt if (vzdt.lt.0) then dQup=Cpfact*(T(i,j,k+1)-T(i,j,k)+dpTdzp) if ((-vzdt).gt.dz) then tmp(i,j,k)=Q(i,j,k)+dQup else tmp(i,j,k)=Q(i,j,k)-dQup*vzdt*dzi endif else dQdwn=Cpfact*(T(i,j,k)-T(i,j,k-1)-dpTdzm) if (vzdt.gt.dz) then tmp(i,j,k)=Q(i,j,k)-dQdwn else tmp(i,j,k)=Q(i,j,k)-dQdwn*vzdt*dzi endif endif 30 continue c We have to do the bottom boundary correctly kmin=kmn(i,j) dpTdzp=-0.5*(rlp(kmin+1)+rlp(kmin))*dz vzdt=vz(i,j,kmin)*dt if (vzdt.lt.0) then dQup=Cpfact*(T(i,j,kmin+1)-T(i,j,kmin)+dpTdzp) if ((-vzdt).gt.dz) then tmp(i,j,kmin)=Q(i,j,kmin)+dQup else tmp(i,j,kmin)=Q(i,j,kmin)-dQup*vzdt*dzi endif else tmp(i,j,kmin)=Q(i,j,kmin) endif 20 continue 10 continue do 40 i=1,nxmx do 50 j=1,nymx kminp=kmn(i,j)+1 do 60 z=kminp,nzmx Q(i,j,k)=tmp(i,j,k) T(i,j,k)=Cpfacti*Q(i,j,k) 60 continue kmin=kmn(i,j) Q(i,j,kmin)=tmp(i,j,kmin) T(i,j,kmin)=Cpfacti*Q(i,j,kmin) 50 continue 40 continue return end c------------------End Qadvectz subroutine -------------- c------------------Begin Vadvectz subroutine -------------- subroutine Vadvectz(v,vz,tmp,nxmx,nymx,nzmx,nxd,nyd,nzd, 1 dz,dt) dimension v(nxd,nyd,nzd),vz(nxd,nyd,nzd) dimension tmp(nxd,nyd,nzd) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) dzi=1/dz do 10 i=1,nxmx do 20 j=1,nymx kminp=kmn(i,j)+1 do 30 k=kminp,nzmx-1 vzdt=vz(i,j,k)*dt if (vzdt.lt.0) then dVup= v(i,j,k+1)-v(i,j,k) if ((-vzdt).gt.dz) then tmp(i,j,k)=v(i,j,k)+dVup else tmp(i,j,k)=v(i,j,k)-dVup*vzdt*dzi endif else dVdwn=v(i,j,k)-v(i,j,k-1) if (vzdt.gt.dz) then tmp(i,j,k)=v(i,j,k)-dVdwn else tmp(i,j,k)=v(i,j,k)-dVdwn*vzdt*dzi endif endif 30 continue 20 continue 10 continue do 40 i=1,nxmx do 50 j=1,nymx kminp=kmn(i,j)+1 do 60 z=kminp,nzmx-1 v(i,j,k)=tmp(i,j,k) 60 continue v(i,j,nzmx)=v(i,j,nzmx-1) 50 continue 40 continue return end c------------------End Vadvectz subroutine -------------- c------------------Begin thadvect------------------------------ subroutine thadvect character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 c Here we had checked the time step, but this is now taken care of c in initval. dtmx=dt dxi=1./dx dyi=1./dy do 20 j=1,nymx do 10 i=2,nxmx-1 vxdt=vx(i,j,nzmx)*dtmx if (vxdt.lt.0) then if ((-vxdt).gt.dx) then tmp(i,j,1)=th(i+1,j) else tmp(i,j,1)=(th(i,j)*(dx+vxdt)-th(i+1,j)*vxdt)*dxi endif else if (vxdt.gt.dx) then tmp(i,j,1)=th(i-1,j) else tmp(i,j,1)=(th(i,j)*(dx-vxdt)+th(i-1,j)*vxdt)*dxi endif endif 10 continue c Now lets do some of the left and right boundaries c first the left boundary i=1 vxdt=vx(i,j,nzmx)*dtmx if (vxdt.lt.0) then if ((-vxdt).gt.dx) then tmp(i,j,1)=th(i+1,j) else tmp(i,j,1)=(th(i,j)*(dx+vxdt)-th(i+1,j)*vxdt)*dxi endif else tmp(i,j,1)=th(i,j) endif c now the the right boundary i=nxmx vxdt=vx(i,j,nzmx)*dtmx if (vxdt.gt.0) then tmp(i,j,1)=th(i,j) else if (vxdt.gt.dx) then tmp(i,j,1)=th(i-1,j) else tmp(i,j,1)=(th(i,j)*(dx-vxdt)+th(i-1,j)*vxdt)*dxi endif endif 20 continue c we load up the new values for the y-advection do 15 i=1,nxmx do 25 j=1,nymx th(i,j)=tmp(i,j,1) 25 continue 15 continue c Now the y-advection loop do 40 i=1,nxmx do 50 j=2,nymx-1 vydt=vy(i,j,nzmx)*dtmx if (vydt.le.0) then if ((-vydt).gt.dy) then tmp(i,j,1)=th(i,j+1) else tmp(i,j,1)=(th(i,j)*(dy+vydt)-th(i,j+1)*vydt)*dyi endif else if (vydt.gt.dy) then tmp(i,j,1)=th(i,j-1) else tmp(i,j,1)=(th(i,j)*(dy-vydt)+th(i,j-1)*vydt)*dyi endif endif 50 continue c Now lets do some of the left and right boundaries c first the bottom boundary j=1 vydt=vy(i,j,nzmx)*dtmx if (vydt.lt.0) then if ((-vydt).gt.dy) then tmp(i,j,1)=th(i,j+1) else tmp(i,j,1)=(th(i,j)*(dy+vydt)-th(i,j+1)*vydt)*dyi endif else tmp(i,j,1)=th(i,j) endif c now the the top boundary j=nymx vydt=vy(i,j,nzmx)*dtmx if (vydt.gt.0) then tmp(i,j,1)=th(i,j) else if (vydt.gt.dy) then tmp(i,j,1)=th(i,j-1) else tmp(i,j,1)=(th(i,j)*(dy-vydt)+th(i,j-1)*vydt)*dyi endif endif c end north and south boundaries 40 continue 60 continue c we load up the new values for the y-advection do 45 i=1,nxmx do 55 j=1,nymx th(i,j)=tmp(i,j,1) 55 continue 45 continue return end c------------------End thadvect------------------------------ c------------------Begin advectxy------------------------------ subroutine advectxy(F,vx,vy,tmp,nxmx,nymx,nzmx,nxd,nyd,nzd, 1 dx,dy,dt,case) dimension F(nxd,nyd,nzd),vx(nxd,nyd,nzd),vy(nxd,nyd,nzd) dimension tmp(nxd,nyd,nzd) common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) c Here we had checked the time step, but this is now taken care of c in initval. dtmx=dt dxi=1./dx dyi=1./dy do 30 k=nzmx-1,2,-1 do 20 j=1,nymx do 10 i=2,nxmx-1 c see if there is solid land on any side: if(k.ge.kmn(i+1,j)) then Fxp=F(i+1,j,k) else Fxp=F(i,j,k)*case endif if(k.ge.kmn(i-1,j)) then Fxm=F(i-1,j,k) else Fxm=F(i,j,k)*case endif c end side land check vxdt=vx(i,j,k)*dtmx if (vxdt.lt.0) then if ((-vxdt).gt.dx) then tmp(i,j,k)=Fxp else tmp(i,j,k)=(F(i,j,k)*(dx+vxdt)-Fxp*vxdt)*dxi endif else if (vxdt.gt.dx) then tmp(i,j,k)=Fxm else tmp(i,j,k)=(F(i,j,k)*(dx-vxdt)+Fxm*vxdt)*dxi endif endif 10 continue c Now lets do some of the left and right boundaries c first the left boundary i=1 vxdt=vx(i,j,k)*dtmx if (vxdt.lt.0) then if ((-vxdt).gt.dx) then tmp(i,j,k)=F(i+1,j,k) else tmp(i,j,k)=(F(i,j,k)*(dx+vxdt)-F(i+1,j,k)*vxdt)*dxi endif else tmp(i,j,k)=F(i,j,k) endif c now the the right boundary i=nxmx vxdt=vx(i,j,k)*dtmx if (vxdt.gt.0) then tmp(i,j,k)=F(i,j,k) else if (vxdt.gt.dx) then tmp(i,j,k)=F(i-1,j,k) else tmp(i,j,k)=(F(i,j,k)*(dx-vxdt)+F(i-1,j,k)*vxdt)*dxi endif endif 20 continue 30 continue c we load up the new values for the y-advection do 15 i=1,nxmx do 25 j=1,nymx kminp=kmn(i,j)+1 do 35 k=kminp,nzmx-1 F(i,j,k)=tmp(i,j,k) 35 continue 25 continue 15 continue c Now the y-advection do 60 k=nzmx-1,2,-1 do 40 i=1,nxmx do 50 j=2,nymx-1 c see if there is solid land on any side: if(k.ge.kmn(i,j+1)) then Fyp=F(i,j+1,k) else Fyp=F(i,j,k)*case endif if(k.ge.kmn(i,j-1)) then Fym=F(i,j-1,k) else Fym=F(i,j,k)*case endif c end side land check vydt=vy(i,j,k)*dtmx if (vydt.le.0) then if ((-vydt).gt.dy) then tmp(i,j,k)=Fyp else tmp(i,j,k)=(F(i,j,k)*(dy+vydt)-Fyp*vydt)*dyi endif else if (vydt.gt.dy) then tmp(i,j,k)=Fym else tmp(i,j,k)=(F(i,j,k)*(dy-vydt)+Fym*vydt)*dyi endif endif 50 continue c Now lets do some of the left and right boundaries c first the bottom boundary j=1 vydt=vy(i,j,k)*dtmx if (vydt.lt.0) then if ((-vydt).gt.dy) then tmp(i,j,k)=F(i,j+1,k) else tmp(i,j,k)=(F(i,j,k)*(dy+vydt)-F(i,j+1,k)*vydt)*dyi endif else tmp(i,j,k)=F(i,j,k) endif c now the the top boundary j=nymx vydt=vy(i,j,k)*dtmx if (vydt.gt.0) then tmp(i,j,k)=F(i,j,k) else if (vydt.gt.dy) then tmp(i,j,k)=F(i,j-1,k) else tmp(i,j,k)=(F(i,j,k)*(dy-vydt)+F(i,j-1,k)*vydt)*dyi endif endif c end north and south boundaries 40 continue 60 continue c we load up the new values for the y-advection do 45 i=1,nxmx do 55 j=1,nymx kminp=kmn(i,j)+1 do 65 k=kminp,nzmx-1 F(i,j,k)=tmp(i,j,k) 65 continue 55 continue 45 continue return end c------------------End subroutine advectxy ----------- c------------------Start subroutine momadvect ----------- subroutine momadvect character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 real tmp3(245,245,95) do 10 i=1,nxmx do 20 j=1,nymx kmin=kmn(i,j) do 30 k=kmin,nzmx tmp2(i,j,k)=vx(i,j,k) tmp3(i,j,k)=vy(i,j,k) 30 continue 20 continue 10 continue case=0. call advectxy(tmp2,vx,vy,tmp,nxmx,nymx,nzmx,245,245,95,dx,dy,dt, 1 case) call Vadvectz(tmp2,vz,tmp,nxmx,nymx,nzmx,245,245,95,dz,dt) case=0. call advectxy(tmp3,vx,vy,tmp,nxmx,nymx,nzmx,245,245,95,dx,dy,dt, 1 case) call Vadvectz(tmp3,vz,tmp,nxmx,nymx,nzmx,245,245,95,dz,dt) do 15 i=1,nxmx do 25 j=1,nymx kmin=kmn(i,j) do 35 k=kmin,nzmx vx(i,j,k)=tmp2(i,j,k) vy(i,j,k)=tmp3(i,j,k) 35 continue if (kmin.ge.2) then do 37 k=1,kmin-1 vy(i,j,k)=0. vx(i,j,k)=0. 37 continue endif 25 continue 15 continue return end c------------------End subroutine momadvect ----------- c------------------Start subroutine accelerate ----------- subroutine accelerate character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 Dnsinv=1./Dens0 do 10 i=1,nxmx do 20 j=1,nymx tmp(i,j,nzmx+1)=0. tmp2(i,j,nzmx+1)=0. kmin=kmn(i,j) do 30 k=nzmx-1,kmin,-1 tmp(i,j,k)= -dt*grdpx(i,j,k) tmp(i,j,nzmx+1)=tmp(i,j,nzmx+1)+tmp(i,j,k) tmp2(i,j,k)= -dt*grdpy(i,j,k) tmp2(i,j,nzmx+1)=tmp2(i,j,nzmx+1)+tmp2(i,j,k) 30 continue tmp(i,j,nzmx+1)=tmp(i,j,nzmx+1)/float(nzmx-kmin+1) tmp2(i,j,nzmx+1)=tmp2(i,j,nzmx+1)/float(nzmx-kmin+1) tmp(i,j,nzmx)=0. tmp2(i,j,nzmx)=0. c kmin=kmn(i,j) kmin=kmn(i,j)+1 c Note that we will need to add a term that helps zero out c the top boundary verticle velocity. do 40 k=nzmx-1,kmin,-1 vxold=vx(i,j,k) vx(i,j,k)=vx(i,j,k)+Dnsinv*tmp(i,j,k) 1 -Dnsinv*tmp(i,j,nzmx+1) 2 +dt*fc(j)*vy(i,j,k) 3 -dt*GPx0*Dnsinv vy(i,j,k)=vy(i,j,k)+Dnsinv*tmp2(i,j,k) 1 -Dnsinv*tmp2(i,j,nzmx+1) 2 -dt*fc(j)*vxold 3 -dt*GPy0*Dnsinv 40 continue c vx(i,j,nzmx)=1.25*vx(i,j,nzmx-1)-0.25*vx(i,j,nzmx-5) c vy(i,j,nzmx)=1.25*vy(i,j,nzmx-1)-0.25*vy(i,j,nzmx-5) vx(i,j,nzmx)=vx(i,j,nzmx-1) vy(i,j,nzmx)=vy(i,j,nzmx-1) 20 continue 10 continue return end c------------------End subroutine accelerate ----------- c------------------Start subroutine heatinputs ----------- subroutine heatinputs character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 pi=4.*atan(1.) sinval=sin((time-6)*pi/12.) Rad=Rad0*(0.5*(sinval+abs(sinval))) c 1 +Rnt do 10 i=1,nxmx do 20 j=1,nymx kmin=kmn(i,j) elev=(kmin)*dz+Z0-dep(i,j) T0=Tz0+dt0(i,j)+realz(elev)*st0 Rad2=Rnt*(T(i,j,kmin)-T0) if(ntmp(i,j).eq.1) then Q(i,j,kmin)=Cpfact*Tcool else if(nsea(i,j).eq.1) then y=float(j-nymx/2)/float(nymx) x=float(i-nxmx/2)/float(nxmx) ymx=0.5*((y-x)-abs(y-x)) Q(i,j,kmin)=Cpfact*(Tsea 1 +dTsea*(((1+(ymx))*(1-(ymx)))-0.5)) c 1 +dTsea*(y-x)) else Q(i,j,kmin)=Q(i,j,kmin) 1 +dt*(Rad*abland-Rad2)*area endif endif T(i,j,kmin)=Q(i,j,kmin)*Cpfacti c Set temperature to specific default if it is underground if(kmin.ge.2) then do 30 k=1,kmin-1 elev=k*dz T(i,j,k)=0-0.01*elev Q(i,j,k)=Cpfact*T(i,j,k) 30 continue endif 20 continue 10 continue return end c------------------End subroutine heatinputs ----------- c------------------Start subroutine momdiffuse ----------- subroutine momdiffuse character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 c First we calculate the stability condition and set our time c step to satisfy it cmax = (2*Cvzmx/(dz*dz)+2*Cvxy1/(dx*dy)+2*Cvxy2/(dx*dy)) dtmx = 0.5/cmax ndt = 1+int(dt/dtmx) dtmx = dt/ndt c Now we loop through with the smaller time step do 5 n=1,ndt cdzmx = dtmx*Cvzmx/(dz*dz) cdzmn = dtmx*Cvzmn/(dz*dz) c The way ground diffusion works, there should be only c one dz... cdzgd = dtmx*Cvzgd/(dz) cdzsea = dtmx*Cvzsea/(dz) cdxy1 = dtmx*Cvxy1/(dx*dy) cdxy2 = dtmx*Cvxy2/(dx*dy) vmix=5. do 10 i=2,nxmx-1 do 20 j=2,nymx-1 kminp=kmn(i,j)+1 do 30 k=kminp,nzmx-1 stratup=(T(i,j,k+1)-T(i,j,k))/dz stratdwn=(T(i,j,k)-T(i,j,k-1))/dz rlpup=0.5*(rlp(k+1)+rlp(k)) rlpdwn=0.5*(rlp(k)+rlp(k-1)) c We need to split the verticle diffusion into an up-component c and a down-component if (stratup.le.rlpup) then cdzup=cdzmx else cdzup=cdzmn endif if (stratdwn.le.rlpdwn) then cdzdn=cdzmx else cdzdn=cdzmn endif c see if there is solid land on any side: if(k.ge.kmn(i,j+1)) then vxyp=vx(i,j+1,k) vyyp=vy(i,j+1,k) else vxyp=vx(i,j,k) vyyp=vy(i,j,k) endif if(k.ge.kmn(i,j-1)) then vxym=vx(i,j-1,k) vyym=vy(i,j-1,k) else vxym=vx(i,j,k) vyym=vy(i,j,k) endif if(k.ge.kmn(i+1,j)) then vxxp=vx(i+1,j,k) vyxp=vy(i+1,j,k) else vxxp=vx(i,j,k) vyxp=vy(i,j,k) endif if(k.ge.kmn(i-1,j)) then vxxm=vx(i-1,j,k) vyxm=vy(i-1,j,k) else vxxm=vx(i,j,k) vyxm=vy(i,j,k) endif c end side land check sumdif=(cdzup+cdzdn+2*cdxy1+2*cdxy2) tmp(i,j,k)=vx(i,j,k)*(1-sumdif) 1 +cdzdn*vx(i,j,k-1)+cdzup*vx(i,j,k+1) 2 +cdxy2*(vxyp+vxym) 3 +cdxy1*(vxxp+vxxm) tmp2(i,j,k)=vy(i,j,k)*(1-sumdif) 1 +cdzdn*vy(i,j,k-1)+cdzup*vy(i,j,k+1) 2 +cdxy1*(vyyp+vyym) 3 +cdxy2*(vyxp+vyxm) 30 continue c Now do the lower boundary.... c Enforce zero velocity at lower boundary kmin=kmn(i,j) c First we use a stability-based criterion stratup=(T(i,j,kmin+1)-T(i,j,kmin))/dz rlpup=0.5*(rlp(kmin+1)+rlp(kmin)) if (stratup.le.rlpup) then cdzup=cdzmx else cdzup=cdzmn endif c Now we will use a criterion based on if it is over the sea if (nsea(i,j).eq.1) then cdzdn=cdzsea else cdzdn=cdzgd endif vamp2=sqrt(vx(i,j,kmin)**2+vy(i,j,kmin)**2) amp=exp(-cdzdn*vamp2) tmp(i,j,kmin)=amp*(vx(i,j,kmin)*(1-cdzup) 1 +cdzup*vx(i,j,kmin+1)) tmp2(i,j,kmin)=amp*(vy(i,j,kmin)*(1-cdzup) 1 +cdzup*vy(i,j,kmin+1)) 20 continue 10 continue c Now we load the new values into the old array... do 40 i=2,nxmx-1 do 50 j=2,nymx-1 kmin=kmn(i,j) do 60 k=kmin,nzmx-1 vx(i,j,k)=tmp(i,j,k) vy(i,j,k)=tmp2(i,j,k) 60 continue if (kmin.ge.2) then do 65 k=1,kmin-1 vx(i,j,k)=0 vy(i,j,k)=0 65 continue endif c vx(i,j,nzmx)=1.25*vx(i,j,nzmx-1)-0.25*vx(i,j,nzmx-5) c vy(i,j,nzmx)=1.25*vy(i,j,nzmx-1)-0.25*vy(i,j,nzmx-5) vx(i,j,nzmx)=vx(i,j,nzmx-1) vy(i,j,nzmx)=vy(i,j,nzmx-1) 50 continue 40 continue c Now we do the side boundaries do 70 k=1,nzmx do 80 i=1,nxmx vx(i,1,k)=vx(i,2,k) vy(i,1,k)=vy(i,2,k) vx(i,nymx,k)=vx(i,nymx-1,k) vy(i,nymx,k)=vy(i,nymx-1,k) 80 continue do 90 j=1,nymx vx(1,j,k)=vx(2,j,k) vy(1,j,k)=vy(2,j,k) vx(nxmx,j,k)=vx(nxmx-1,j,k) vy(nxmx,j,k)=vy(nxmx-1,j,k) 90 continue 70 continue 5 continue return end c------------------End subroutine momdiffuse ----------- c------------------Start subroutine heatdiffuse ----------- subroutine heatdiffuse character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 c First we calculate the stability condition and set our time c step to satisfy it cmax = (2*Czmx/(dz*dz)+4*Cxy/(dx*dy)) dtmx = 0.5/cmax ndt = 1+int(dt/dtmx) dtmx = dt/ndt c Now we loop through with the smaller time step do 5 n=1,ndt c Now begin the heat diffusion loop cdzmx = dtmx*Czmx/(dz*dz) cdzmn = dtmx*Czmn/(dz*dz) cdzgd = dtmx*Czgd/(dz*dz) cdxy = dtmx*Cxy/(dx*dy) vmix=5. do 10 i=2,nxmx-1 do 20 j=2,nymx-1 c Q(i,j,nzmx)=1.25*Q(i,j,nzmx-1)-0.25*Q(i,j,nzmx-5) c T(i,j,nzmx)=1.25*T(i,j,nzmx-1)-0.25*T(i,j,nzmx-5) c T(i,j,nzmx)=T(i,j,nzmx-1) c Q(i,j,nzmx)=Q(i,j,nzmx-1) kminp=kmn(i,j)+1 do 30 k=kminp,nzmx-1 dpQp=-dz*Cpfact*0.5*(rlp(k+1)+rlp(k)) dpQm=-dz*Cpfact*0.5*(rlp(k)+rlp(k-1)) c see if there is solid land on any side: if(k.ge.kmn(i,j+1)) then Qyp=Q(i,j+1,k) else Qyp=Q(i,j,k) endif if(k.ge.kmn(i,j-1)) then Qym=Q(i,j-1,k) else Qym=Q(i,j,k) endif if(k.ge.kmn(i+1,j)) then Qxp=Q(i+1,j,k) else Qxp=Q(i,j,k) endif if(k.ge.kmn(i-1,j)) then Qxm=Q(i-1,j,k) else Qxm=Q(i,j,k) endif c end side land check stratup=(T(i,j,k+1)-T(i,j,k))/dz stratdwn=(T(i,j,k)-T(i,j,k-1))/dz rlpup=0.5*(rlp(k+1)+rlp(k)) rlpdwn=0.5*(rlp(k)+rlp(k-1)) if (stratup.le.rlpup) then cdzup=cdzmx else cdzup=cdzmn endif if (stratdwn.le.rlpdwn) then cdzdn=cdzmx else cdzdn=cdzmn endif sumdif=(cdzup+cdzdn+4*cdxy) tmp(i,j,k)=Q(i,j,k)*(1-sumdif) 1 +cdzdn*(Q(i,j,k-1)-dpQm)+cdzup*(Q(i,j,k+1)+dpQp) 2 +cdxy*(Qxp+Qxm 3 +Qyp+Qym) 30 continue c Now do the lower boundary.... kmin=kmn(i,j) dpQp=-dz*Cpfact*0.5*(rlp(kmin+1)+rlp(kmin)) stratup=(T(i,j,kmin+1)-T(i,j,kmin))/dz rlpup=0.5*(rlp(kmin+1)+rlp(kmin)) if (stratup.le.rlpup) then cdz=cdzmx else cdz=cdzmn endif sumdif=cdz tmp(i,j,kmin)=Q(i,j,kmin)*(1-sumdif) 1 +cdz*(Q(i,j,kmin+1)+dpQp) Tg(i,j)=T(i,j,kmin) 1 +Cpfacti*(tmp(i,j,kmin)-Q(i,j,kmin))/(dtmx*dzgd*Czgd) 20 continue 10 continue c Now we load the new values into the old array... do 40 i=2,nxmx-1 do 50 j=2,nymx-1 kmin=kmn(i,j) do 60 k=kmin,nzmx-1 Q(i,j,k)=tmp(i,j,k) T(i,j,k)=Q(i,j,k)*Cpfacti 60 continue c Enforce temperature condition on the sea kmin=kmn(i,j) if (nsea(i,j).eq.1) then y=float(j-nymx/2)/float(nymx) x=float(i-nxmx/2)/float(nxmx) ymx=0.5*((y-x)-abs(y-x)) Q(i,j,kmin)=Cpfact*(Tsea 1 +dTsea*(((1+(ymx))*(1-(ymx)))-0.5)) c 1 +dTsea*(y-x)) T(i,j,kmin)=Q(i,j,kmin)*Cpfacti endif if (ntmp(i,j).eq.1) then T(i,j,kmin)=Tcool Q(i,j,kmin)=Cpfact*Tcool endif c Now lets reinforce the top boundary condition: c Q(i,j,nzmx)=1.25*Q(i,j,nzmx-1)-0.25*Q(i,j,nzmx-5) c T(i,j,nzmx)=1.25*T(i,j,nzmx-1)-0.25*T(i,j,nzmx-5) T(i,j,nzmx)=Tup Q(i,j,nzmx)=Cpfact*T(i,j,nzmx) 50 continue 40 continue c Now we do the side boundaries delT1=-5. delT2=-5. elsea=1500. c We impose a particular marine air profile in two marine basins c top marine basin is between the 208 and 771 kilometer point, c in the x-direction, bottom marine basin is above 1028 kilometer c point in x-direction, and below the 450 kilometer point in the y c direction nxmin1=208000/dx nxmax1=771000/dx nxmin2=963000/dx nymax2=450000/dy do 70 k=1,nzmx elev=(k)*dz+Z0-dep(i,j) elf=(elsea-elev)/elsea do 80 i=1,nxmx j=1 kmin=kmn(i,j) c if ((i.ge.nxmin2).and.(elev.lt.elsea) c 1 .and.(k.gt.kmin)) then c T(i,j,k)=Tz0-delT2*elf+(Tup-Tz0)*((dz*float(k-1)+Z0)/Zup) c Q(i,j,k)=Cpfact*T(i,j,k) c else Q(i,1,k)=Q(i,2,k) T(i,1,k)=T(i,2,k) c endif j=nymx kmin=kmn(i,j) c if ((i.ge.nxmin1).and.(elev.lt.elsea) c 1 .and.(k.gt.kmin).and.(i.le.nxmax1)) then c T(i,j,k)=Tz0-delT1*elf+(Tup-Tz0)*((dz*float(k-1)+Z0)/Zup) c Q(i,j,k)=Cpfact*T(i,j,k) c else Q(i,nymx,k)=Q(i,nymx-1,k) T(i,nymx,k)=T(i,nymx-1,k) c endif 80 continue do 90 j=1,nymx i=1 kmin=kmn(i,j) Q(1,j,k)=Q(2,j,k) T(1,j,k)=T(2,j,k) i=nxmx kmin=kmn(i,j) c if ((j.le.nymax2).and.(elev.lt.elsea) c 1 .and.(k.gt.kmin)) then c T(i,j,k)=Tz0-delT2*elf+(Tup-Tz0)*((dz*float(k-1)+Z0)/Zup) c Q(i,j,k)=Cpfact*T(i,j,k) c else Q(nxmx,j,k)=Q(nxmx-1,j,k) T(nxmx,j,k)=T(nxmx-1,j,k) c endif 90 continue 70 continue 5 continue return end c------------------End subroutine heatdiffuse ----------- c------------------Begin subroutine refine ----------- c This subroutine does the job of refining or interpolating c a field on a grid in a way that minimized curvature, but c retains the average values of the original coarse field. Subroutine refine(rf,x,nxmx,nymx,ndx1,ndy1,ndx2,ndy2,nref) dimension rf(ndx1,ndy1),x(ndx2,ndy2) dimension tmp(490,490) c write(*,*) " In refine " do 10 i=1,nxmx do 20 j=1,nymx do 30 ni=1,nref do 32 nj=1,nref in=nref*(i-1)+ni jn=nref*(j-1)+nj rf(in,jn)=x(i,j) 32 continue 30 continue 20 continue 10 continue c So we have loaded up the average pixel values. Now we c need to smooth and readjust First smooth: nxmx2=nref*nxmx nymx2=nref*nymx nit=2 do 99 it=1,nit do 40 i=2,nxmx2-1 do 50 j=2,nymx2-1 tmp(i,j)=0.25*(rf(i+1,j)+rf(i-1,j)+rf(i,j+1)+rf(i,j-1)) 50 continue tmp(i,1)=0.25*(rf(i+1,1)+rf(i-1,1)+rf(i,2)+rf(i,1)) j=nymx2 tmp(i,j)=0.25*(rf(i+1,j)+rf(i-1,j)+rf(i,j-1)+rf(i,j)) 40 continue do 41 j=2,nymx2-1 tmp(1,j)=0.25*(rf(1,j+1)+rf(1,j-1)+rf(2,j)+rf(1,j)) i=nxmx2 tmp(i,j)=0.25*(rf(i,j+1)+rf(i,j-1)+rf(i-1,j)+rf(i,j)) 41 continue tmp(1,1)=rf(1,1) tmp(1,nymx2)=rf(1,nymx2) tmp(nxmx2,1)=rf(nxmx2,1) tmp(nxmx2,nymx2)=rf(nxmx2,nymx2) do 60 i=1,nxmx2 do 70 j=1,nymx2 rf(i,j)=tmp(i,j) 70 continue 60 continue if (it.ne.nit) then c Now lets do the average pixel adjustment step nav=nref**2 do 110 i=1,nxmx do 120 j=1,nymx avrf=0 do 130 ni=1,nref do 132 nj=1,nref in=nref*(i-1)+ni jn=nref*(j-1)+nj avrf=avrf+rf(in,jn) 132 continue 130 continue avrf=avrf/nav dref=x(i,j)-avrf do 135 ni=1,nref do 137 nj=1,nref in=nref*(i-1)+ni jn=nref*(j-1)+nj rf(in,jn)=rf(in,jn)-dref 137 continue 135 continue 120 continue 110 continue endif 99 continue return end c------------------End subroutine refine ----------- c------------------Begin subroutine initcase ----------- c This subroutine initializes the parameter values for the c different cases being run. So that we can run lots of cases, c and see what parameter values give reasonable answers subroutine initcase(nc) character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 c Reset time time=6. c Set some of the parameter values: if (nc.eq.1) then c GPx0=0.000 c GPy0=0.000 dTsea=14. dTatm=6. Tz0=14. Tsea=27.-dTsea*0.5 st0=-0.0060 outname="dataout1.txt" endif if (nc.eq.2) then close(unit=11) c GPx0=0.0001 c GPy0=0.0001 dTsea=8. dTatm=6. Tz0=18. Tsea=28.-dTsea*0.5 st0=-0.0075 outname="dataout2.txt" endif if (nc.eq.3) then close(unit=11) c GPx0=0.000 c GPy0=0.00015 dTsea=-4. dTatm=6. Tz0=25. Tsea=30.-dTsea*0.5 st0=-0.0090 outname="dataout3.txt" endif if (nc.eq.4) then close(unit=11) c GPx0=-0.0001 c GPy0=0.0001 dTsea=18. dTatm=6. Tz0=18. Tsea=28.-dTsea*0.5 st0=-0.0080 outname="dataout4.txt" endif if (nc.eq.5) then close(unit=11) GPx0=-0.00015 GPy0=0.000 outname="dataout5.txt" endif if (nc.eq.6) then close(unit=11) GPx0=-0.0001 GPy0=-0.0001 outname="dataout6.txt" endif if (nc.eq.7) then close(unit=11) GPx0=0.000 GPy0=-0.00015 outname="dataout7.txt" endif if (nc.eq.8) then close(unit=11) GPx0=0.0001 GPy0=-0.0001 outname="dataout8.txt" endif if (nc.eq.9) then close(unit=11) GPx0=0.00015 GPy0=0.000 outname="dataout9.txt" endif c Now open the datafile for output open(unit=11,file=outname) return end c------------------End subroutine initcase ----------- c------------------Start subroutine radloss ----------- subroutine radloss character*12 outname common/arrays/ Q(245,245,95),T(245,245,95), 1 vx(245,245,95),vy(245,245,95),vz(245,245,95), 2 tmp(245,245,95),tmp2(245,245,95), 3 grdpx(245,245,95),grdpy(245,245,95) c 4 ,wc(245,245,95) common/runpar/itmx,nxmx,nymx,nzmx,ndit,dt,dx,dy,dz,nx0,ny0, 1 nz0,nxsea,area,outname common/kmin/kmn(245,245),nsea(245,245),ntmp(245,245) common/elev/slpx(245,245),slpy(245,245),dep(245,245),dt0(245,245) 1 ,Tg(245,245),th(245,245) common/physpar/Cd,Dens0,Rad0,Rnt,Tup,Tsea,Tz0,Z0,Zup,rlp(95),st0, 1 time,Czmx,Czmn,Cxy,Cvzmx,Cvzmn,Czgd,Cvzsea,Cvzgd,Cp,Cpsea, 2 abland,absea,Cpfact,Cpfacti,alpha,fc(245),Tcool,GPx0,GPy0,dTsea 3 ,zmax,dzgd,dTatm,Cvxy1,Cvxy2 c We expect a general radiative loss from the atmosphere, and a c daytime radiative gain. The delta T for the equilibrium temperature c profile is set in initval tloss=30*3600 tloss2=10*3600 ftloss=exp(-dt/tloss) ftloss2=exp(-dt/tloss2) pi=4.*atan(1.) c dTz0 is the amplitude of the drive, Tresp is the amplitude c of the equilibrium response of bulk radiative forcing. c To enforce the daily bulk temperature oscillation, we calculate c the derivative of the temperature varition and add that time c the time step to the existing temperature. We set the temperature c oscillation to peak at 4 pm and nadir at 4 am. Tresp=1.5 c Derivative: freq=2.*pi/24 dTresp=Tresp*freq*cos(freq*(time-10))/(3600) c Note that we also have to dissipate bulk momentum in order c to prevent a linear net momentum production instability. c we will set bulk momentum dissipation to the same time scale do 10 i=1,nxmx do 20 j=1,nymx kmin=kmn(i,j) do 30 k=kmin+1,nzmx elev=(k)*dz+Z0 T0=Tz0+dTatm+realz(elev)*st0 c 1 +dt0(i,j) T(i,j,k)=T0+dTresp*dt+(T(i,j,k)-T0)*ftloss Q(i,j,k)=Cpfact*T(i,j,k) vx(i,j,k)=vx(i,j,k)*ftloss vy(i,j,k)=vy(i,j,k)*ftloss c We also want to dampen boundary velocities to decrease inflow c instabilities if ((i.le.2).or.(i.ge.nxmx-1).or. 1 (j.le.2).or.(j.ge.nymx-1)) then vx(i,j,k)=vx(i,j,k)*ftloss2 vy(i,j,k)=vy(i,j,k)*ftloss2 c T(i,j,k)=T0+(T(i,j,k)-T0)*ftloss2 c Q(i,j,k)=Cpfact*T(i,j,k) endif 30 continue c For the surface deviation, we want to also dissipate waves on that c surface, otherwise it will collect energy without dissipation th(i,j)=ftloss2*th(i,j) 20 continue 10 continue return end c------------------End subroutine radloss -----------