I receive the following error
Compiling file: tropic.f
Warning: Extension: Tab character in format at (1)
C:\Users\Marchant\Desktop\tropic.f(432) : error - Expected a right parenthesis in expression at column 72
Warning: Rank mismatch in argument 'tk' at (1) (scalar and rank-1)
Warning: Rank mismatch in argument 't' at (1) (scalar and rank-1)
Warning: Rank mismatch in argument 'tk' at (1) (scalar and rank-1)
Warning: Rank mismatch in argument 't' at (1) (scalar and rank-1)
Compilation failed.
in this program,
dimension ts1(3),ts2(3),ta1(3),ta2(3),out(14,300)
real lwc, lambda
common /params/xkap,pr,tr,xl,cp,rgas,grav,taue,taus,taup,tauc
common /param2/ pbot,ptop,dp,gam, bt,ct,tao,a21,lambda,lwc
common /heat/ beta,olr1,olr2,alb0,albgr,expo1,expo2,alb1,alb2
pbot=1.0e5
ptop=2.0e4
dp=pbot-ptop
open(12,file='tropic.in',form='formatted')
read(12,*) itermx, delt, iprint
read(12,*) lambda, gam, bt, ct, a1
read(12,*) beta,olr1,olr2,alb0,albgr,expo1,expo2
write(*,*) 'olr1=',olr1,', olr2=',olr2,', expo1=',expo1,', expo2='
1 ,expo2
c ** Set relative areas of convecting a1 and nonconvecting a2 regions.
c a1=.3
tao=265.
alpha=0.06
alpha2=alpha/2.
alpha1=1.-alpha
c expo1=80.
c expo2=80.
expa1=0.
expa2=0.
co=4.2e7
ca=1.0e7
xkap=0.288
rvap=461.
cp=1004.
rgas=287.
grav=9.81
c gam=1.0e-3
c lambda=1.0e3
pr=1.0e5
tr=300.
xl=2.5e6
write(*,*) ' gam=',gam
c** structure of output array
c out(1)=a1; 2=gam; 3=lambda
c 4=ts1 5=ts2 6=alb1 7=alb2
c 8=r1 9=r2 10=ts1tend 11=ts2tend
c 13=thet1 14=thet2
ikase=0
c ********* BIG LOOP ****************
do 888 nn=1,2
a1=0.1+0.2*nn
do 888 ll=1,7
c gam=1.0e-3*facg
gam=1/1024.*2.0**(ll-1)
do 888 mm=1,7
c lambda=1.0e+3*facl
lambda=64*2.0**(mm-1)
c write(*,*) '*******************************'
c write(*,*) 'GAM=',gam,', LAMBDA=',lambda,', A1=',a1
a2=1.-a1
a21=a2/a1
a12=a1/a2
c initialize variables
do i = 1,3
ts1(i)=301.
ts2(i)=300.
ta1(i)=302.
ta2(i)=300.
end do
is=1
js=2
tdelto=2.*delt/co
tdelta=2.*delt/ca
c write(*,999) ts1(js),ts2(js),ta1(js),ta2(js),r1,r2,ra1,ra2
999 format(1x,9f8.1)
c write(*,*) pbot,ptop,dp,pr,gam,bt,ct,tao,a21,lambda,lwc
ikase=ikase+1
c*** Time Loop *****
do 1000 it=1,itermx
dta=ta1(js)-ta2(js)
dto=ts1(js)-ts2(js)
call radiat(ts1(js),ts2(js),ta1(js),ta2(js),r1,r2,ra1,ra2)
call theta(ts1(js),ts2(js),ta1(js),ta2(js),demdp,demd2,deddp)
c** Note that demdp = del(theta)/grav
ts1(3)=ts1(is)+tdelto*(r1-gam*dto*cp*demdp-expo1)
ts2(3)=ts2(is)+tdelto*(r2+a12*gam*dto*cp*demdp-expo2)
c ta1(3)=ta1(is)+tdelta*(ra1-a21*gam*dto*cp*demdp-expa1)
c ta2(3)=ta2(is)+tdelta*(ra2+gam*dto*cp*deddp-expa2)
c apply Robert/Asselin filter
ts1(js)=ts1(js)*alpha1 +alpha2*(ts1(3)+ts1(is))
ts2(js)=ts2(js)*alpha1 +alpha2*(ts2(3)+ts2(is))
c if((it-1)/iprint*iprint.eq.it-1) then
if((it.eq.itermx)) then
time=(it-1)*delt/86400.
ts1tend=(r1-gam*dto*cp*demdp-expo1)*86400./co
ts2tend=(r2+a12*gam*dto*cp*demdp-expo2)*86400./co
c ta1tend=(-a21*gam*dto*cp*demdp)
c ta2tend=( gam*dto*cp*demdp)
thet1=thet(ts1,qsat(ts1,pbot),pbot)
thet2=thet(ts2,qsat(ts2,pbot),pbot)
c** structure of output array
c out(1)=a1; 2=gam; 3=lambda
c 4=ts1 5=ts2 6=alb1 7=alb2
c 8=r1 9=r2 10=ts1tend 11=ts2tend
c 12=thet1 13=thet2
c Set up array
out(1,ikase)=a1
out(2,ikase)=gam
out(3,ikase)=lambda
out(4,ikase)=ts1(js)
out(5,ikase)=ts2(js)
out(6,ikase)=alb1
out(7,ikase)=alb2
out(8,ikase)=r1
out(9,ikase)=r2
out(10,ikase)=ts1tend
out(11,ikase)=ts2tend
out(12,ikase)=thet1
out(13,ikase)=thet2
out(14,ikase)=qsat(ts1(js),pr)
c write(*,*) 'Day=',time, ', iter=',it
c write(*,*) a21,gam,dto,cp,demdp
c write(*,*) 'demdp, demd2,deddp', demdp, demd2,deddp
c write(*,*) 'lwc=',lwc,alb1, alb2
c*********x*********x*********x*********x*********x*********x*********x**********
c write(*,*) ' ts1, ts2, ta1, ta2, r1, r2, ra1,
c 1 ra2'
c write(*,999) ts1(3),ts2(3),ta1(3),ta2(3),r1,r2,ra1,ra2
c write(*,999) ts1(js),ts2(js),ta1(js),ta2(js),r1,r2,ra1,ra2
c write(*,998) ts1tend,ts2tend,ta1tend,ta2tend, thet1, thet2
998 format(1x,8f10.5)
endif
c ** Update Variables
is=3-is
js=3-js
ts1(js)=ts1(3)
ts2(js)=ts2(3)
ta1(js)=ta1(3)
ta2(js)=ta2(3)
1000 continue
888 continue
open(13,file='tropic.out',form='formatted')
c*********x*********x*********x*********x*********x*********x*********x**********
write(*,*) ' A1 gam lambda ts1 ts2 alb1
1alb2 r1 r2 ts1tend ts2tend thet1 thet2 qsat'
write(13,*) ' A1 gam lambda ts1 ts2 alb1
1alb2 r1 r2 ts1tend ts2tend thet1 thet2 qsat'
do ii=1,ikase
xkrap=out(2,ii)*out(3,ii)
write(*,789) (out(j,ii),j=1,14),xkrap
write(13,789) (out(j,ii),j=1,14),xkrap
789 format(1x,f6.1,f9.5,7f9.2,2f9.5,2f8.2,2f8.4)
enddo
stop
end
c ******************************************************
subroutine theta(ts1,ts2,ta1,ta2,demdp,demd2,deddp)
c ** This subroutine finds the theta gradients
real lwc, lambda
common /param2/ pbot,ptop,dp,gam, bt,ct,tao,a21,lambda,lwc
common /params/xkap,pr,tr,xl,cp,rgas,grav,taue,taus,taup,tauc
demdp=(thet(ts1,qsat(ts1,pbot),pbot)-thet(ts2,qsat(ts2,pbot),
1 pbot))/9.81
c 1 pbot))/dp
demd2=(thet(ta1,0.001,ptop)-thet(ts1,qsat(ts1,pbot),pbot))
1 /9.81
c 1 /dp
deddp=(thet(ts1,0.00001,ptop)-thet(ts2,0.00001,pbot))/9.81
c 1 /dp
return
end
c ******************************************************
subroutine radiat(ts1,ts2,ta1,ta2,r1,r2,ra1,ra2)
real lwc, lambda
common /param2/ pbot,ptop,dp,gam, bt,ct,tao,a21,lambda,lwc
common /params/xkap,pr,tr,xl,cp,rgas,grav,taue,taus,taup,tauc
common /heat/ beta,olr1,olr2,alb0,albgr,expo1,expo2,alb1,alb2
dta=ta1-ta2
dto=ts1-ts2
if(dto.gt.0.0) then
c ** radiation parameterization for atmosphere
ra1=-40-bt*(ta1-tao)+ct*(ts1-(ta1+29))
ra2=-200-bt*(ta2-tao)+ct*(ts2-(ta2+29))
c ** Get liquid water content
c lwc=lambda*a21*gam*abs(dto)*qsat(ts1,pr)
c ** Get albedo as function of LWC
alb2=alb0
alb1=alb0+lambda*gam*abs(dto)*qsat(ts1,pr)
if(alb1.gt.0.75) alb1=0.75
r1=400.*(1.-alb1)-olr1-beta*(ts1-300.)
r2=400.*(1.-alb2)-olr2-beta*(ts2-300.)
else
c ** here ts2 is hotter than ts1
c ** radiation parameterization for atmosphere
ra1=-200-bt*(ta1-tao)+ct*(ts1-(ta1+29))
ra2=-40-bt*(ta2-tao)+ct*(ts2-(ta2+29))
c ** Get liquid water content
c lwc=lambda*gam*abs(dto)*qsat(ts2,pr)
c ** Get albedo as function of LWC
alb1=alb0
alb2=alb0+lambda*gam*abs(dto)*qsat(ts2,pr)
if(alb2.gt.0.75) alb2=0.75
r1=400.*(1.-alb1)-olr2-beta*(ts1-300.)
r2=400.*(1.-alb2)-olr1-beta*(ts2-300.)
endif
c write(*,*) 'lwc=',lwc,', alb1,2=',alb1,alb2,', r,ra-',r1,r2,ra1,ra2
return
end
c*********x*********x*********x*********x*********x*********x*********x**********
c*************************************************************
function temp(the,rv,p)
c** Function calculates temperature given thetaE, rv and p
common /params/xkap,pr,tr,xl,cp,rgas,grav,taue,taus,taup,tauc
temp=the/((pr/p)**xkap*exp(xl*rv/(cp*tr)))
return
end
c*************************************************************
function thet(t,rv,p)
c** Function calculates thetaE given t, rv and p
common /params/xkap,pr,tr,xl,cp,rgas,grav,taue,taus,taup,tauc
thet=t*(pr/p)**xkap*exp(xl*rv/(cp*tr))
return
end
c*************************************************************
function thets(t,p)
c** Function calculates thetaEsaturate given t and p
common /params/xkap,pr,tr,xl,cp,rgas,grav,taue,taus,taup,tauc
if(t.lt.273.15) then
es=esice(t)
else
es=esat(t)
endif
rs=0.622*es/(p-es)
thets=t*(pr/p)**xkap*exp(xl*rs/(cp*tr))
return
end
c*************************************************************
subroutine plevs(p,xlp,dlp,dp)
c** Subroutine to set pressure levels
parameter(ilx=25)
dimension p(ilx),xlp(ilx),dlp(ilx),dp(ilx)
write(*,*) 'Setting Pressure Levels'
write(*,*) ' i p(i) dp(i) logp dlogp'
pmin=2000.
pmax=101300.
delpo=pmax-pmin
delp=delpo/(ilx-1)
do i=1,ilx
p(i)=pmin+(i-1.)*delp
xlp(i)=alog(p(i))
end do
do i=1,ilx-1
dlp(i)=xlp(i+1)-xlp(i)
dp(i)=p(i+1)-p(i)
end do
dlp(ilx)=0.0
do i=1,ilx
write(*,*) i,p(i),dp(i),xlp(i),dlp(i)
end do
return
end
c*************************************************************
subroutine radini(teq,p,t,sst)
c** Calculates variables needed by radiation relaxation code
parameter (ilx=25)
dimension p(ilx),t(ilx),teq(ilx)
do i=1,ilx
if(p(i).lt.12000.) then
teq(i)=t(i)
c elseif(p(i).gt.80000.) then
else
teq(i)=t(i)-10.
c teq(i)=t(i)-(p(ilx)/10000.)*2.
endif
end do
return
end
c*************************************************************
subroutine initlz(the,rt,rs,t,rv,p,sst)
c** Subroutine to set initial values of all variables
parameter (ilx=25)
dimension the(ilx),rt(ilx),rs(ilx),t(ilx),rv(ilx),
1 p(ilx)
common /params/xkap,pr,tr,xl,cp,rgas,grav,taue,taus,taup,tauc
ttrop=200.
tsurf=300.
ptrop=10000.
dtdp=(tsurf-ttrop)/(p(ilx)-ptrop)
relhum=0.80
c** Set T(p)
do i=1,ilx
if(p(i).lt.ptrop) then
t(i)=200.+10.*(ptrop-p(i))/(ptrop-p(1))
else
t(i)=200.+dtdp*(p(i)-ptrop)
endif
end do
c** Next calculate vapor mixing ratio and thetaE
write(*,*) 'index, pressure, temp., vapor mr, thetaE'
do i=1,ilx
if(p(i).lt.ptrop) then
rfrac=0.05
else
rfrac=relhum
endif
if(t(i).lt.273.) then
es=esice(t(i))
else
es=esat(t(i))
endif
rv(i)=rfrac*0.622*es/(p(i)-es)
rs(i)=0.622*es/(p(i)-es)
rt(i)=rv(i)
the(i)=t(i)*(pr/p(i))**xkap*exp(xl*rv(i)/(cp*tr))
write(*,100) i,p(i),t(i),rv(i),the(i)
100 format(1x,i3,f12.1,f7.1,e13.3,f7.1)
end do
return
end
c*************************************************************
function signum(x)
c** Hankel function
if(x.eq.0) then
signum=1.
else
signum=(abs(x)+x)*0.5/abs(x)
endif
return
end
c*************************************************************
subroutine zero(x,n)
dimension x(n)
do i=1,n
x(i)=0.0
end do
return
end
C#######################################################################
FUNCTION ESICE(TK)
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE WITH RESPECT TO
C ICE ESICE (Pascals) GIVEN THE TEMPERATURE T (Kelvin). DLH 11.19.97
C THE FORMULA USED IS BASED UPON THE INTEGRATION OF THE CLAUSIUS-
C CLAPEYRON EQUATION BY GOFF AND GRATCH. THE FORMULA APPEARS ON P.350
C OF THE SMITHSONIAN METEOROLOGICAL TABLES, SIXTH REVISED EDITION,
C 1963.
DATA CTA,EIS/273.15,6.1071/
C CTA = DIFFERENCE BETWEEN KELVIN AND CELSIUS TEMPERATURE
C EIS = SATURATION VAPOR PRESSURE (MB) OVER A WATER-ICE MIXTURE AT 0C
DATA C1,C2,C3/9.09718,3.56654,0.876793/
C C1,C2,C3 = EMPIRICAL COEFFICIENTS IN THE GOFF-GRATCH FORMULA
c**** Convert to Celsius
c tc=t-273.15
IF (TK.LE.CTA) GO TO 5
ESICE = 99999.
WRITE(6,3)ESICE
3 FORMAT(' SATURATION VAPOR PRESSURE FOR ICE CANNOT BE COMPUTED',
1 /' FOR TEMPERATURE > 0C. ESICE =',F7.0)
RETURN
5 CONTINUE
C FREEZING POINT OF WATER (K)
TF = CTA
C GOFF-GRATCH FORMULA
RHS = -C1*(TF/TK-1.)-C2*ALOG10(TF/TK)+C3*(1.-TK/TF)+ALOG10(EIS)
ESI = 10.**RHS
IF (ESI.LT.0.) ESI = 0.
ESICE = ESI*100.
RETURN
END
C#######################################################################
FUNCTION ESAT(TK)
C THIS FUNCTION RETURNS THE SATURATION VAPOR PRESSURE OVER
C WATER (Pa) GIVEN THE TEMPERATURE (Kelvin). DLH 11.19.97
C THE ALGORITHM IS DUE TO NORDQUIST, W.S.,1973: "NUMERICAL APPROXIMA-
C TIONS OF SELECTED METEORLOLGICAL PARAMETERS FOR CLOUD PHYSICS PROB-
C LEMS," ECOM-5475, ATMOSPHERIC SCIENCES LABORATORY, U.S. ARMY
C ELECTRONICS COMMAND, WHITE SANDS MISSILE RANGE, NEW MEXICO 88002.
IF (TD.NE. 99999.0) THEN
C IF (TD.NE.-1001.0) THEN
c**** Convert to Celsius
c TK = TD+273.15
P1 = 11.344-0.0303998*TK
P2 = 3.49149-1302.8844/TK
C1 = 23.832241-5.02808*ALOG10(TK)
ESAT = 100.*10.**(C1-1.3816E-7*10.**P1+8.1328E-3*10.**P2-2949.076/TK)
else
esat = 0.
END IF
RETURN
END
C#######################################################################
function qsat(tk,p)
qsat=esat(tk)*0.622/p
return
end
Can someone show me how to fix this? its a fortran77 file being compiled in mingw gfortran
