PROGRAM TWOD C****************************************************************************** C This program runs a two-dimensional hydrostatic model containing C parameterizations of convection, clouds and radiative transfer. The model C may be oriented in the east-west direction with periodic boundaries, C or in the north-south direction with rigid boundaries. Radiation and C ocean mixed layer temperature may be calculated or specified in the input C sounding file. The model is nonlinear, except that the advection of the C overturning circulation by itself is neglected. Standard Fickian C diffusion modifies the wind variables, but not temperature or specific C humidity. The parameter NLO specifies the number of grid columns; it C must be at least 4. The number of vertical levels is specified in the C input sounding file. C C There are two input files that must be supplied: sounding.in, which C contains one vertical sounding of all the model variables for each C grid column, and params.in, which contains most of the parameters C used to control the model. There must be at least as many sounding C columns in sounding.in as NLO-2. C C Output consists of the file sounding.out, whose format is identical to C that of sounding.in, and a series of matrix files containing time- C distance series and distance-pressure cross-section averaged over C the last n time steps of the integration, where n is specified in C params.in. One may specify these files to be in matlab-readable ASCII C format or comma-separated ASCII format. The matlab script files C imenux and imenuy provide ways of displaying the output. C C The main program and subroutine CONVECT were created by Kerry Emanuel (Emanuel C and Zivkovic-Rothman, J. Atmos. Sci., 1999, 1766-1782), the fractioanl cloud C scheme by Sandrine Bony (Bony and Emanuel, J. Atmos. Sci., 2001, 3158-3183) C and the radiation code was developed by J.-J. Morcrette (J. Geophys. Res., 1991, C 9121-9132). The radiation subroutine may be deleted if radiative cooling is C specified; this will result in faster compilation and running. C C***************************************************************************** C C Note: The parameter NA must be greater than or equal to the total number C of model levels at which the state variables are defined. This will C equal the number of tropospheric levels plus 9. NOTE THAT NA MUST BE C SET TO THE SAME VALUE IN SUBROUTINE CONVECT. NTIME is at least the C number of time steps that certain quanitities like precipitation are C output to the graphics arrays. NCL is the number of cloud types C used by the radiation code. C PARAMETER (NA=50, NAP=NA+1, NTIME=5000, NCL=1) PARAMETER (NLO=4) C C Dimension various arrays C C Temperatures: C REAL T1(NLO,NA), T2(NLO,NA), T3(NLO,NA), TTEMP(NA) REAL TCONV(NA),RCONV(NA), T0(NLO,NA) REAL TVP(NA), TMEAN(NLO,NA) C C Surface temperatures: C REAL TS1(NLO),TS2(NLO),TS3(NLO) C C Specific humidities: C REAL R1(NLO,NA), R2(NLO,NA), R3(NLO,NA), RS(NA) REAL RT(NLO,NA), RV2(NA) REAL RTEMP(NA), R0(NLO,NA), RMEAN(NLO,NA) C C Other quantities in convection scheme: C REAL UC(NA),VC(NA),VCO(NA),UCO(NA) REAL TRA1(NA,1) C C Pressures: C REAL P(NA), PH(NAP), PHRAD(NAP) REAL SPONGE(NA) C C Latitudes: C REAL RLAT(NLO), COSZA(NLO), COSLAT(NLO) C C Ozone used in radiative calculation: C REAL O3(NA) C C Streamfunction and zonal wind C REAL PSI(NLO,NA), U(NA), FC(NLO) REAL ETA1(NLO,NA),ETA2(NLO,NA),ETA3(NLO,NA) REAL U1(NLO,NA), U2(NLO,NA), U3(NLO,NA) C C Background meridional wind and T gradient C REAL VZ(NA), TX(NA) C C *** Applied Omega *** C REAL AOMEGA(NA) C C Forcings: C REAL FT(NA), FR(NA), FTNET(NLO,NA), FRNET(NLO,NA) REAL FTS(NLO),ETAF(NLO,NA),UF(NLO,NA),RADFLUX(NLO) REAL FTRA(NA,1),FUC(NA),FVC(NA) REAL FTADJ(NA), FRADJ(NA),VMAXK(NLO) C C Radiative arrays: C REAL RCOOL(NLO,NA), CS(NCL), PCTOP(NCL), PCBOT(NCL), TAUS(NCL) REAL FNETIR(NAP), FNETSO(NAP), COOLR(NA), HEATR(NA), HOUR REAL QCONDT(NA) REAL HOUR0 REAL MEAN,QSAT,CLDF(NA),CLDQ(NA),SIGSUB(NA),QSUBGRID : ,gsol,minimum,xo,ho,pi,qmin,qmax,pas,signe,x : ,delta(NA),TIMEL,essai : ,RSNEWLS(NA), RNEWLS(NA), TNEWLS(NA) : ,CLDFRAD(NA), CLDQRAD(NA) REAL rmumean(NLO), fracmean(NLO), albmean(NLO) CHARACTER*1 RADINT, TSINT, TDEP, DDEP, DARAD, ANRAD, CALB CHARACTER*1 RESTART C C Clouds and radiation C integer julien real CLDEMI(NA), CLDTAU(NA), CLDFICE(NA), CLDT, CLDWP : ,zlongi, dist, RLON, gmtime, albsol, RCOTN(NLO,NA) : ,pplay(NA), paprs(NAP), RCOTNOLD(NA) : ,RCLW(NA),RCSW(NA),RCLWOLD(NA),RCSWOLD(NA) : ,fract,rmu0 : ,heat(NA),heat0(NA),cool(NA),cool0(NA) : ,radsol,albpla,topsw,toplw,solsw,sollw : ,topsw0,toplw0,solsw0,sollw0 : ,radsolold,radsoln C C Surface arrays C REAL EFRAC(NLO), DEPTH(NLO), SALB(NLO), FSOIL(NLO) C C Soil arrays C REAL SM1(NLO), SM2(NLO), SM3(NLO) C C Convective momentum transport C REAL FVCNET(NLO,NA) C C Saved arrays C REAL RV2M(NA),CLDEMIM(NA),CLDTAUM(NA),CLDFRADM(NA) C C Character constants C CHARACTER*1 WDEP,BCS,FWRITE,ICLOUD,IVAPOR,ISOIL C C Convective quantities used for diagnostics: C REAL MP(NA), RP(NA), SIJ(NA,NA), WATER(NA) REAL CBMFA(NLO),FRSUB(NA) REAL MUP(NA), MDOWN(NA), ENT(NA), DET(NA) REAL FRDET(NA) C C Time-dependent quantities for graphics output: C REAL QGRAPH(NLO,NTIME),OGRAPH(NLO,NTIME),PGRAPH(NLO,NTIME) REAL VGRAPH(NLO,NTIME),UTGRAPH(NLO,NTIME),GTIME(NTIME) REAL SSTGRAPH(NLO,NTIME), MPIGRAPH(NLO,NTIME) REAL T100GRAPH(NLO,NTIME), T100BAR(NLO), EFBAR(NLO) REAL SWGRAPH(NLO,NTIME), CAPEGRAPH(NLO,NTIME) REAL CAPEBAR(NLO), QRBAR(NLO), QRGRAPH(NLO,NTIME) REAL OLRGRAPH(NLO,NTIME),HCGRAPH(NLO,NTIME) REAL LCGRAPH(NLO,NTIME), EFRACGRAPH(NLO,NTIME) C C Space-dependent quantities for output: C REAL OLGRAPH(NLO,NA),MUPGRAPH(NLO,NA),MDGRAPH(NLO,NA) REAL PREBAR(NLO),RHGRAPH(NLO,NA),ULGRAPH(NLO,NA) REAL OBAR(NLO),VBAR(NLO),UBAR(NLO),SSTBAR(NLO) REAL PSIGRAPH(NLO,NA),MUPG(NLO,NA),MDG(NLO,NA),RHG(NLO,NA) REAL OLG(NLO,NA),TEMPGRAPH(NLO,NA), AMGRAPH(NLO,NA) REAL AFLUXGRAPH(NLO),OFLUXGRAPH(NLO),AFLUX(NLO),OFLUX(NLO) REAL VOCEAN(NLO), CLDTEMP(NLO,NA), CLDGRAPH(NLO,NA) REAL NETFLUX(NLO), NETFLUXGRAPH(NLO) REAL SWBAR(NLO), TEMPGRAPHNET(NLO,NA) REAL OLRBAR(NLO),LCBAR(NLO),HCBAR(NLO) C C Special graphics arrays: C REAL DISX(NLO) C C Latent heat at T=0 C C REAL LV0 C LOGICAL ok, correc, correc2, rad1st, radcomp : ,newobs,nodomega,newqlw,evap_prec : ,correctd C C Common block with subroutine CONVECT: C COMMON / CVT / TVP,MP,RP,SIJ,MUP,MDOWN,ENT,DET,FRDET,FRSUB, 1 WATER,INB, CAPEM C C Common block with subroutine CLOUDS_SUB_LS: C REAL CLDL(NA),CLDA(NA),CLDK(NA),CLDS(NA) COMMON / cldpara / CLDL, CLDA, CLDK, CLDS c ------------------------------------------------------- c OPTIONS OF THE SIMULATION: c ------------------------------------------------------- c if radcomp=T: the rad tendency added to net tendencies is computed: radcomp = .TRUE. c maximum precipitation efficiency: EPMAX = 0.999 c minimum precipitation efficiency: EPMIN = 0.0 C rad1st = .FALSE. C C Open input and output files: C OPEN(UNIT=8, FILE='params.in', STATUS='OLD') OPEN(UNIT=11, FILE='output/error.out', STATUS='unknown') C C C Read in the parameters of the model from the file params.in C READ(8,*) READ(8,*) READ(8,*) READ(8,*) READ(8,*)RESTART READ(8,*)ENDTIME READ(8,*)DT READ(8,*) READ(8,*) READ(8,*) READ(8,*)AWN READ(8,*)RLATB READ(8,*)RLATE READ(8,*)RLONB READ(8,*)RLONE READ(8,*)KSOIL READ(8,*)BCS READ(8,*) READ(8,*) READ(8,*) READ(8,*)RADINT READ(8,*)RADFREQ READ(8,*)TSINT READ(8,*)IVAPOR READ(8,*)ICLOUD READ(8,*)SCON READ(8,*)CO2 READ(8,*)MONTH READ(8,*)IDAY READ(8,*)HOUR READ(8,*)TDEP READ(8,*)DDEP READ(8,*)DARAD READ(8,*)ANRAD READ(8,*)CALB READ(8,*)ALB READ(8,*)ALBLAND READ(8,*) READ(8,*) READ(8,*) READ(8,*)WDEP READ(8,*)US0 READ(8,*)VS0 READ(8,*)WS0 READ(8,*) READ(8,*) READ(8,*) READ(8,*)ISOIL READ(8,*)EFRACR READ(8,*)TAUSOIL READ(8,*)TAURUNOFF READ(8,*)SBETA READ(8,*) READ(8,*) READ(8,*) READ(8,*)DEPTHR READ(8,*)CCOUP READ(8,*) READ(8,*) READ(8,*) READ(8,*)DPBL READ(8,*)DAMP READ(8,*)SCOEF READ(8,*)PFIX READ(8,*) READ(8,*) READ(8,*) READ(8,*)PRINFREQ READ(8,*)PAVE READ(8,*)VTRANS READ(8,*)FWRITE C C Renormalize certain parameters C C RDAMP=4./(RADFREQ*3600.) AOMEGMAX=0.0 RDAMP=0.0 DT=DT*60.0 DAMP=1.0/(24.*3600.*DAMP) PNUM=100.*DAMP SCOEF=1.0/(60.*SCOEF) TDEEP=TDEEP+273.15 AOMEGMAX=AOMEGMAX/36.0 ALENGTH=1.0E6*AWN DX=ALENGTH/FLOAT(NLO-2) DXI=1.0/DX DTI=1.0/DT DPBLI=1.0/DPBL ENDTIME=ENDTIME*3600.0*24.0+HOUR*3600. RADFREQ=RADFREQ*3600.0 PRINFREQ=PRINFREQ*3600.0 PAVE=PAVE*24.*3600. PAVE=MIN(PAVE,ENDTIME) PAFRAC=DT/(PAVE+1.0) WS0=WS0*WS0 IF(RESTART.EQ.'Y')RESTART='y' IF(RADINT.EQ.'Y')RADINT='y' IF(RADINT.EQ.'N')RADINT='n' IF(RADINT.EQ.'n')radcomp=.FALSE. IF(TSINT.EQ.'Y')TSINT='y' IF(TSINT.EQ.'N')TSINT='n' IF(ANRAD.EQ.'Y')ANRAD='y' IF(ANRAD.EQ.'N')ANRAD='n' IF(TDEP.EQ.'Y')TDEP='y' IF(TDEP.EQ.'N')TDEP='n' IF(DDEP.EQ.'Y')DDEP='y' IF(DDEP.EQ.'N')DDEP='n' IF(CALB.EQ.'Y')CALB='y' IF(CALB.EQ.'N')CALB='n' IF(ANRAD.EQ.'y')THEN TDEP='n' DDEP='n' END IF IF(DARAD.EQ.'y')THEN TDEP='n' END IF MONTH0=MONTH IDAY0=IDAY HOUR0=HOUR WFAC=0.0 KQMAX=0 IF(MONTH.EQ.1)THEN JULIENDAY=IDAY ELSE IF(MONTH.EQ.2)THEN JULIENDAY=31+IDAY ELSE IF(MONTH.EQ.3)THEN JULIENDAY=31+28+IDAY ELSE IF(MONTH.EQ.4)THEN JULIENDAY=31+28+31+IDAY ELSE IF(MONTH.EQ.5)THEN JULIENDAY=31+28+31+30+IDAY ELSE IF(MONTH.EQ.6)THEN JULIENDAY=31+28+31+30+31+IDAY ELSE IF(MONTH.EQ.7)THEN JULIENDAY=31+28+31+30+31+30+IDAY ELSE IF(MONTH.EQ.8)THEN JULIENDAY=31+28+31+30+31+30+31+IDAY ELSE IF(MONTH.EQ.9)THEN JULIENDAY=31+28+31+30+31+30+31+31 +IDAY ELSE IF(MONTH.EQ.10)THEN JULIENDAY=31+28+31+30+31+30+31+31+30+IDAY ELSE IF(MONTH.EQ.11)THEN JULIENDAY=31+28+31+30+31+30+31+31+30+31+IDAY ELSE IF(MONTH.EQ.12)THEN JULIENDAY=31+28+31+30+31+30+31+31+30+31+30+IDAY END IF C C Latitude and Corilois parameter C PIM=ACOS(-1.0)/180. DO K=1,NLO RLAT(K)=RLATB+(RLATE-RLATB)*FLOAT(K-1)/FLOAT(NLO-1) FC(K)=1.454E-4*SIN(PIM*RLAT(K)) COSLAT(K)=COS(PIM*RLAT(K)) END DO C C Calculate annual average earth-sun distance, solar zenith angle and albedo C IF(ANRAD.EQ.'y')THEN distm=0.0 do k=1,nlo rmumean(k)=0.0 fracmean(k)=0.0 albmean(k)=0.0 end do DO JDAY=1,365 CALL orbite(FLOAT(JDAY),zlongi,dist) distm=distm+dist do k=2,nlo rlatx=0.5*(RLAT(K)+RLAT(K-1)) CALL angle(zlongi,rlatx,fract,rmu0) CALL alboc(FLOAT(JDAY),rlatx,albsol) rmumean(k)=rmumean(k)+rmu0 fracmean(k)=fracmean(k)+fract albmean(k)=albmean(k)+albsol end do END DO distm=distm/365. do k=2,nlo rmumean(k)=rmumean(k)/365. fracmean(k)=fracmean(k)/365. albmean(k)=albmean(k)/365. end do END IF C C Ocean mixed layer depth and effective soil depth C DO K=1,NLO IF(K.LT.KSOIL)THEN DEPTH(K)=DEPTHR ELSE DEPTH(K)=1.0 END IF END DO C C Soil moisture parameters C TAUSOIL=TAUSOIL*24.*3600. TAURUNOFF=TAURUNOFF*24.*3600. C TAUSOILI=1./TAUSOIL TAURUNOFFI=1./TAURUNOFF IF(TAURUNOFF.GT.100.0)THEN TAURUNOFFI=0.0 END IF C C Read in the initial soundings from the file sounding.in C IF(RESTART.EQ.'y')THEN OPEN(UNIT=9, FILE='output/sounding.out', STATUS='OLD') ELSE OPEN(UNIT=9, FILE='sounding.in', STATUS='OLD') END IF C DO 42 J=2,NLO-1 READ(9,20)NP,TS1(J) 20 FORMAT(T15,I3,T55,F5.2,////) c c TS1(J)=27.5-2.5*TANH(4.*(float(J-10)/FLOAT(NLO-3))) c TS1(J)=TS1(J)+273.15 TS2(J)=TS1(J) DO 40 I=1,NP READ(9,30)P(I),T1(J,I),R1(J,I),O3(I),PSI(J,I),RCOOL(J,I), 1 U1(J,I) 30 FORMAT(3X,F6.1,T12,F8.3,T25,F10.6,T39,F6.3,T47,F6.2, 1 T58,F7.3,T64,F5.1) C C Renormalize certain quantities and produce reverse order arrays for C radiative subroutine C c PSI(J,I)=0.0 c U1(J,I)=0.0 C T1(J,I)=T1(J,I)+273.15 R1(J,I)=0.001*R1(J,I) T2(J,I)=T1(J,I) T0(J,I)=T1(J,I) TMEAN(J,I)=T1(J,I) R2(J,I)=R1(J,I) R0(J,I)=R1(J,I) RMEAN(J,I)=R1(J,I) RCOOL(J,I)=RCOOL(J,I)/(24.*3600.) c PSI(J,I)=0.0 PSI(J,I)=1.0E4*PSI(J,I) IF(P(I).LT.(PFIX+1.0))THEN PSI(J,I+1)=0.0 U1(J,I)=0.0 END IF U2(J,I)=U1(J,I) C IF(J.EQ.(NLO-1))THEN T1(J,I)=T1(J,I)+0.1 T2(J,I)=T1(J,I) END IF C 40 CONTINUE 42 CONTINUE C CLOSE(UNIT=8) CLOSE(UNIT=9) C C Background zonal wind C DO I=1,NP c VZ(I)=20.0*(P(1)-P(I))/(P(1)-500.0) VZ(I)=0.0 VZ(I)=MIN(VZ(I),20.0) c vz(i)=vz(i)+5.0 END DO DO I=2,NP-1 c TX(I)=P(I)*FC(2)*(VZ(I+1)-VZ(I-1))/(287.*(P(I-1)- c 1 P(I+1))) TX(I)=0.0 END DO TX(1)=TX(2) TX(NP)=TX(NP-1) C c TS1(2)=TS1(2)+0./(1.+0.15*9.0) c TS2(2)=TS1(2) C C Surface pressure and surface exchange coefficient C PG=1012.0 CD=1.5E-3 C C Asselin Filter Value C AFILT=0.2 C C Renormalize ozone, define levels C PT=200.0 NPTOP=1 N=1 DO 43 I=1,NP O3(I)=1.0E-6*O3(I) IF(P(I).GE.(PFIX-1.0))NPTOP=MAX(NPTOP,I) IF(P(I).GE.78.0)N=MAX(N,I) IF(I.GT.1)THEN PH(I)=0.5*(P(I)+P(I-1)) PHRAD(NP+2-I)=PH(I) END IF UC(I)=0.0 VC(I)=0.0 TRA1(I,1)=0.0 C C Applied Omega C AOMEGA(I)=4.*AOMEGMAX*(P(1)-P(I))*(P(I)-PT)/ 1 ((P(1)-PT)**2) IF(P(I).LT.PT)AOMEGA(I)=0.0 C C Sponge layer coefficients C SPONGE(I)=SCOEF*(1.-P(I)/85.0) SPONGE(I)=MAX(SPONGE(I),0.0) C 43 CONTINUE C C NPTOP=MIN(NPTOP,NP-1) N=MIN(N,(NP-1)) C PH(1)=2.*P(1)-PH(2) PHRAD(NP+1)=PH(1) PH(NP+1)=2.*PH(NP)-PH(NP-1) PH(NP+1)=MAX(0.1,PH(NP+1)) PHRAD(1)=PH(NP+1) C DPH1I=1.0/(PH(1)-PH(2)) DPHNPI=1.0/(PH(NP)-PH(NP+1)) DP1I=1.0/(P(1)-P(2)) DPNPI=1.0/(P(NP-1)-P(NP)) C C Assign thermodynamic constants C CPD=1005.7 CPV=1870.0 CL=2500.0 RD=287.04 RV=461.5 LV0=2.501E6 G=9.8 ROWL=1000.0 EXPAN=3.2E-4 EPS=RD/RV EPSI=1./EPS CPVMCL=CL-CPV C C Initialize certain quantities and arrays C QSPEED=0.0 FS=0.0 FL=0.0 DO 45 J=2,NLO-1 FTS(J)=0.0 PSI(J,1)=0.0 ETA1(J,1)=0.0 ETA2(J,1)=0.0 PSI(J,NP)=0.0 ETA1(J,NP)=0.0 ETA2(J,NP)=0.0 RADFLUX(J)=0.0 CBMFA(J)=0.0 C SM1(J)=0.0 SM2(J)=SM1(J) SALB(J)=ALB 45 CONTINUE C DO 60 I=1,NP COOLR(I)=0.0 HEATR(I)=0.0 QCONDT(I)=0.0 CLDF(I)=0.0 CLDQ(I)=0.0 SIGSUB(I)=0.0 RV2M(I)=0.0 CLDFRADM(I)=0.0 CLDEMIM(I)=0.0 CLDTAUM(I)=0.0 DO 51 K=2,NLO-1 UF(K,I)=0.0 IF(I.NE.1.AND.I.NE.NP)THEN ETA1(K,I)=-1.0E-4*((PSI(K,I+1)-PSI(K,I))/(PH(I+1)-PH(I))- 1 (PSI(K,I)-PSI(K,I-1))/(PH(I)-PH(I-1)))/(P(I)-P(I-1)) ETA2(K,I)=ETA1(K,I) END IF ETAF(K,I)=0.0 FTNET(K,I)=0.0 FRNET(K,I)=0.0 FVCNET(K,I)=0.0 MUPGRAPH(K,I)=0.0 MDGRAPH(K,I)=0.0 OLGRAPH(K,I)=0.0 RHGRAPH(K,I)=0.0 ULGRAPH(K,I)=0.0 PSIGRAPH(K,I)=0.0 TEMPGRAPH(K,I)=0.0 AMGRAPH(K,I)=0.0 CAPEGRAPH(K,I)=0.0 CLDTEMP(K,I)=0.0 CLDGRAPH(K,I)=0.0 51 CONTINUE 60 CONTINUE DO 65 K=2,NLO-1 PREBAR(K)=0.0 VBAR(K)=0.0 UBAR(K)=0.0 SSTBAR(K)=0.0 T100BAR(K)=0.0 SWBAR(K)=0.0 CAPEBAR(K)=0.0 QRBAR(K)=0.0 OLRBAR(K)=0.0 LCBAR(K)=0.0 HCBAR(K)=0.0 EFBAR(K)=0.0 65 CONTINUE DO K=1,NLO AFLUXGRAPH(K)=0.0 OFLUXGRAPH(K)=0.0 NETFLUXGRAPH(K)=0.0 AFLUX(K)=0.0 OFLUX(K)=0.0 VOCEAN(K)=0.0 NETFLUX(K)=0.0 VMAXK(K)=0.0 END DO IF((ICLOUD.EQ.'n'.OR.IVAPOR.EQ.'n').AND.FWRITE.NE.'y')THEN OPEN(10,FILE='output/variables',form='unformatted', 1 status='unknown') READ(10) RV2M,CLDFRADM,CLDEMIM,CLDTAUM CLOSE(10) END IF C C Initialize counters C TIME=HOUR*3600. RADTIME=0.0 PRINTIME=0.0 NT=0 NPR=1 NPRINT=0 TPERR=0.0 OLRERR=0. ASRERR=0. OLRBIA=0. ASRBIA=0. OLRNB=0. ASRNB=0. NRHI=0 NTPI=0 RHERRI=0. TPERRI=0. NSAVE=0 NPRIRAD=0 C C Begin time loop at next statement C 100 CONTINUE C TIME=TIME+DT RADTIME=RADTIME+DT PRINTIME=PRINTIME+DT NT=NT+1 IF(TIME.GT.ENDTIME)GOTO 200 NPRINT=NPRINT+1 c IF(TIME.GT.(50.*24.*3600.))CD=0.0 C c Relate evaporative fraction to soil moisture c DO I=1,NLO-1 IF(I.LT.KSOIL)THEN EFRAC(I)=1.0 ELSE IF(ISOIL.EQ.'y')THEN EFRAC(I)=0.5*(1.+TANH(SBETA*SM2(I))) ELSE EFRAC(I)=EFRACR END IF END IF END DO C C *** Boundary values *** C IF(BCS.EQ.'p')THEN c TS1(1)=TS1(NLO-1) TS2(1)=TS2(NLO-1) TS1(NLO)=TS1(2) TS2(NLO)=TS2(2) SM1(1)=SM1(NLO-1) SM2(1)=SM2(NLO-1) SM1(NLO)=SM1(1) SM2(NLO)=SM2(2) DO 105 I=1,NP U1(1,I)=U1(NLO-1,I) U2(1,I)=U2(NLO-1,I) U1(NLO,I)=U1(2,I) U2(NLO,I)=U2(2,I) T1(1,I)=T1(NLO-1,I) T2(1,I)=T2(NLO-1,I) T1(NLO,I)=T1(2,I) T2(NLO,I)=T2(2,I) R1(1,I)=R1(NLO-1,I) R2(1,I)=R2(NLO-1,I) R1(NLO,I)=R1(2,I) R2(NLO,I)=R2(2,I) c c if(p(i).gt.620.0)then c r1(1,i)=0.004 c r2(1,i)=0.004 c else if(p(i).gt.240.and.p(i).le.620.)then c r1(1,i)=0.0002*(p(i)-200.)/620. c r2(1,i)=0.0002*(p(i)-200.)/620. c else c r1(1,i)=2.0e-6 c r2(1,i)=2.0e-6 c end if c r1(nlo,i)=r1(nlo-1,i) c r2(nlo,i)=r2(nlo-1,i) c PSI(1,I)=PSI(NLO-1,I) PSI(NLO,I)=PSI(2,I) ETA1(1,I)=ETA1(NLO-1,I) ETA1(NLO,I)=ETA1(2,I) FVCNET(1,I)=FVCNET(NLO-1,I) FVCNET(NLO,I)=FVCNET(2,I) 105 CONTINUE C ELSE C TS1(1)=TS1(2) TS2(1)=TS2(2) TS1(NLO)=TS1(NLO-1) TS2(NLO)=TS2(NLO-1) SM1(1)=SM1(2) SM2(1)=SM2(2) SM1(NLO)=SM1(NLO-1) SM2(NLO)=SM2(NLO-1) DO 106 I=1,NP U1(1,I)=0.0 U2(1,I)=0.0 U1(NLO,I)=0.0 U2(NLO,I)=0.0 T1(1,I)=T1(2,I) T2(1,I)=T2(2,I) T1(NLO,I)=T1(NLO-1,I) T2(NLO,I)=T2(NLO-1,I) R1(1,I)=R1(2,I) R2(1,I)=R2(2,I) R1(NLO,I)=R1(NLO-1,I) R2(NLO,I)=R2(NLO-1,I) PSI(1,I)=0.0 PSI(NLO,I)=0.0 ETA1(1,I)=0.0 ETA1(NLO,I)=0.0 PSI(NLO-1,I)=0.0 ETA1(NLO-1,I)=0.0 FVCNET(1,I)=FVCNET(2,I) FVCNET(NLO,I)=FVCNET(NLO-1,I) 106 CONTINUE C END IF SUMFLUX=0.0 DELTGLOBAL=(ABS(TS2(1)-TS2(NLO)))**(1./3) DELTGLOBAL=MAX(DELTGLOBAL,1.0) C DO 155 K=2,NLO-1 C COSLATBAR=0.5*(COSLAT(K)+COSLAT(K-1)) COSLATBARI=1./COSLATBAR FCBAR=0.5*(FC(K)+FC(K-1)) FTS(K)=0.0 IF(TDEP.EQ.'y')THEN XT=TIME/(3600.*24.)+HOUR0/24. HOUR=24.*(XT-AINT(XT)) END IF IF(DDEP.EQ.'y'.OR.DDEP.EQ.'Y')THEN XT=TIME/(3600.*24.)+FLOAT(IDAY0)+30.*FLOAT(MONTH0-1)+ 1 HOUR0/24. DAY=XT-365.*AINT(XT/365.) AMON=1.1+AINT((DAY-1.)/30.) AMON=MIN(12.1,AMON) MONTH=AMON MONTH=MAX(MONTH,1) IDAY=DAY-30.*FLOAT(MONTH-1) END IF C C Calculate saturation specific humidity and some temporary arrays C TC=TS1(K)-273.15 ES=6.112*EXP(17.67*TC/(243.5+TC)) RSS=0.98*EPS*ES/(PG-(1.-EPS)*ES) PRADJ=0.0 TG=TS1(K) RG=RSS C DO 110 I=1,NP FTNET(K,I)=0.0 FRNET(K,I)=0.0 TC=T2(K,I)-273.15 IF(TC.GE.0.0)THEN ES=6.112*EXP(17.67*TC/(243.5+TC)) ELSE ES=EXP(23.33086-6111.72784/T2(K,I)+0.15215*LOG(T2(K,I))) END IF ES=MIN(ES,P(I)) RS(I)=EPS*ES/(P(I)-(1.-EPS)*ES) TTEMP(I)=T2(K,I) RTEMP(I)=R2(K,I) TCONV(I)=T2(K,I) RCONV(I)=R2(K,I) UC(I)=U2(K,I) C 110 CONTINUE C DO I=1,NP-1 VC(I)=0.005*(PSI(K-1,I+1)-PSI(K-1,I)+PSI(K,I+1)-PSI(K,I))/ 1 (PH(I)-PH(I+1))+VZ(I) VCO(I)=VC(I) UCO(I)=UC(I) END DO C CBMF=CBMFA(K) C C Call convection subroutine C CALL CONVECT(TCONV,RCONV,RS,UC,VC,TRA1,P,PH,NA,N,0,DT,IFLAG,FT,FR, 1 FUC,FVC,FTRA,PRECIP,WD,TPRIME,QPRIME,QCONDT,CBMF,FS,FL) C C If necessary, write to error file C IF(IFLAG.EQ.4)THEN WRITE(11,120)NT 120 FORMAT(5X,'CFL condition violated at time step ',I8) END IF C C Store convective momentum flux for vorticity calculation C DO I=1,N FVCNET(K,I)=FVC(I)+(VC(I)-VCO(I))/(0.5*3600.) FUC(I)=FUC(I)+(UC(I)-UCO(I))/(0.5*3600.) END DO C PREBAR(K)=PREBAR(K)+PRECIP CAPEBAR(K)=CAPEBAR(K)+CAPEM QRBAR(K)=QRBAR(K)+RCONV(1) OBAR(K)=OBAR(K)+36.*(COSLAT(K)*PSI(K,9)- 1 COSLAT(K-1)*PSI(K-1,9))*DXI*COSLATBARI VBAR(K)=VBAR(K)+0.01*PSI(K,2)*DPH1I UBAR(K)=UBAR(K)+U2(K,1) SSTBAR(K)=SSTBAR(K)+TS2(K)-273.15 EFBAR(K)=EFBAR(K)+EFRAC(K) T100BAR(K)=T100BAR(K)+T2(K,37)-273.15 C C Add convective tendencies to net tendencies C IF(K.EQ.2)INBA=INB RINB=MIN(R1(1,INBA),8.0E-5) DO 125 I=1,N FTNET(K,I)=FTNET(K,I)+FT(I) FRNET(K,I)=FRNET(K,I)+FR(I) R2(K,I)=RCONV(I) T2(K,I)=TCONV(I) 125 CONTINUE CBMFA(K)=CBMF C C *** Relax water vapor above tropopause to tropopause value *** C DO 126 I=INB+1,NP FRNET(K,I)=FRNET(K,I)+2.0E-6*(0.4*RS(I)*P(I)/P(INB+1)-R1(K,I)) 126 CONTINUE C C -- sb: C ------------------------------------------------------------------- C C Compute the cloud fraction and the cloud water content: C C 1) make supersaturation adjustment C (precipitate only a fraction of the supersaturation) C C 2) compute the cloud fraction using the in-cloud water C content derived from the convection scheme (the subgrid-scale C part of the in-cloud water) and that derived from the C large-scale condensation scheme (or supersaturation sdjustment). C C ------------------------------------------------------------------- C CALL CLOUDS_SUB_LS_40(NP,RCONV,RS,TCONV,P,PH,DT,qcondt : ,CLDF,CLDQ,PRADJ,FTADJ,FRADJ,SIGSUB) C C Add tendencies due to supersaturation adjustment: C DO 111 I = 1, NP C FRNET(K,I)=FRNET(K,I)+FRADJ(I) FTNET(K,I)=FTNET(K,I)+FTADJ(I) C 111 CONTINUE C C Add precip from supersaturation adjustment to convective precip: C PREBAR(K)=PREBAR(K)+PRADJ C C ------------------------------------------------------------------- C C Compute cloud optical properties: C C ------------------------------------------------------------------- c pressure: mb -> Pa: DO I = 1, NP paprs(I) = PH(I)*100. pplay(I) = P(I)*100. ENDDO paprs(NP+1) = PH(NP+1)*100. c IF(ICLOUD.EQ.'n')THEN DO I=1,NP CLDF(I)=CLDFRADM(I) END DO END IF c AHCSUM=1.0 ALCSUM=1.0 DO I = 1, NP CLDFRAD(I) = CLDF(I) CLDTEMP(K,I)=CLDF(I) CLDQRAD(I) = CLDQ(I) IF(P(I).GT.522.0)THEN ALCSUM=ALCSUM*(1.-CLDF(I)) ELSE AHCSUM=AHCSUM*(1.-CLDF(I)) END IF ENDDO LCBAR(K)=LCBAR(K)+1.-ALCSUM HCBAR(K)=HCBAR(K)+1.-AHCSUM C CALL OPTICAL(NP,TCONV,paprs,pplay,CLDFRAD,CLDQRAD : ,CLDEMI,CLDTAU,CLDFICE,CLDT,CLDWP) C ------------------------------------------------------------------- C C Prepare and call radiation: C C ------------------------------------------------------------------- C IF(RADINT.EQ.'y')THEN IF(TIME.LT.(1.5*DT).OR.RADTIME.GT.(RADFREQ-10.0))THEN IF(K.EQ.(NLO-1))RADTIME=0.0 c c -- local hour (RLON=longitude): c RLON=RLONB+(RLONE-RLONB)*FLOAT(K-1)/FLOAT(NLO-1) TIMEL = TIME + RLON/360.*24.*3600. C C -- Use annual means, if specified: C IF(ANRAD.EQ.'y')THEN dist=distm rmu0=rmumean(k) fract=fracmean(k) albsol=albmean(k) ELSE c c -- Earth-Sun distance: c IF(DDEP.EQ.'y')THEN IF(TDEP.EQ.'n')THEN julien = JULIENDAY + INT(TIMEL)/(3600*24) rjulien=float(julien) ELSE rjulien=float(julienday)+timel/(3600.*24.) END IF ELSEIF(TDEP.EQ.'y')THEN rjulien=float(julienday)+timel/(3600.*24.)- 1 INT(timel/(3600.*24.)) ELSE julien = JULIENDAY rjulien=float(julien) END IF c CALL orbite(rjulien,zlongi,dist) c -- solar zenith angle and surface albedo: c RLATX=0.5*(RLAT(K)+RLAT(K-1)) IF(DARAD.EQ.'y')THEN ! no diurnal cycle CALL angle(zlongi,RLATX,fract,rmu0) c print*, rlatx, fract CALL alboc(rjulien,RLATX,albsol) ELSE IF(TDEP.EQ.'y')THEN gmtime = MOD(TIME,3600.*24.)/86400. - 0.5*RADFREQ/86400. CALL zenang(zlongi,gmtime,RADFREQ,RLATX,RLON,rmu0,fract) ELSE gmtime=HOUR/24.-0.05 CALL zenang(zlongi,gmtime,0.1,RLATX,RLON,rmu0,fract) END IF CALL alboc_cd(rmu0,albsol) ENDIF C END IF c c -- call radiation: IF(CALB.EQ.'n'.AND.K.LT.KSOIL)albsol=ALB IF(K.GE.KSOIL)albsol=ALBLAND DO I = 1, NP RV2(I) = RCONV(I) - CLDQ(I)*CLDF(I) RV2(I) = MAX( MIN(RV2(I),RCONV(I)) , 0.0 ) TMEAN(K,I)=T2(K,I) RMEAN(K,I)=R2(K,I) ENDDO C IF(FWRITE.EQ.'y')THEN NSAVE=NSAVE+1 DO I=1,NP RV2M(I)=RV2M(I)+RV2(I) CLDFRADM(I)=CLDFRADM(I)+CLDFRAD(I) CLDEMIM(I)=CLDEMIM(I)+CLDEMI(I) CLDTAUM(I)=CLDTAUM(I)+CLDTAU(I) END DO END IF IF(ICLOUD.EQ.'n')THEN DO I=1,NP CLDFRAD(I)=CLDFRADM(I) CLDEMI(I)=CLDEMIM(I) CLDTAU(I)=CLDTAUM(I) END DO END IF IF(IVAPOR.EQ.'n')THEN DO I=1,NP RV2(I)=MIN(RV2M(I),RS(I)) END DO END IF C TG=TS2(K) C CALL radlwsw e (dist, rmu0, fract, CO2, SCON, e paprs, pplay,TG,albsol,TCONV,RV2,O3, e CLDFRAD, CLDEMI, CLDTAU, s heat,heat0,cool,cool0,radsol,albpla, s topsw,toplw,solsw,sollw, s topsw0,toplw0,solsw0,sollw0) C C Calculate net radiative cooling rate C SWBAR(K)=SWBAR(K)+solsw OLRBAR(K)=OLRBAR(K)+toplw IF(K.EQ.2)NPRIRAD=NPRIRAD+1 RADFLUX(K)=radsol NETFLUX(K)=NETFLUX(K-1)+COSLATBAR*DX*(topsw-toplw) radsoln=radsol DO 130 I=1,NP RCOTN(K,I)=(cool(I)-heat(I))/(24.*3600.) RCLW(I)=cool0(I)-cool(I) ! diagnostic only RCSW(I)=heat(I)-heat0(I) ! diagnostic only c thermo constants used in Morcrette's radiation are different from those c used in the calling program, so we renormalize the cooling rates: RCOTN(K,I)=RCOTN(K,I)*1004.709/(CPD*(1.-RCONV(I))+ 1 CPV*RCONV(I)) RCOTN(K,I)=RCOTN(K,I)*G/9.80665 130 CONTINUE c if (radcomp) rad1st = .TRUE. c END IF END IF C sb -- C C Add radiative cooling to net heating rate C DO 135 I=1,NP if (RADINT.EQ.'y') then FTNET(K,I)=FTNET(K,I)-RCOTN(K,I) RCOOL(K,I)=RCOTN(K,I) FTNET(K,I)=FTNET(K,I)-RDAMP*(T1(K,I)-TMEAN(K,I)) FRNET(K,I)=FRNET(K,I)-RDAMP*(R1(K,I)-RMEAN(K,I)) else IF(P(I).GT.95.0)THEN c -- sb: cc FTNET(I)=FTNET(I)-RCOT cc RCOOL(I)=RCOT c FTNET(I)=FTNET(I)-RCOTU(I) FTNET(K,I)=FTNET(K,I)-RCOOL(K,I) c RCOOL(I)=RCOTU(I) c sb -- ELSE RCOOL(K,I)=0.0 END IF endif ! radcomp 135 CONTINUE C C Calculate surface flux forcing C TSA=T1(K,1)*(PG/P(1))**(RD/CPD) ROWS=PG/(RD*TSA*(1.+R2(K,1)*(EPSI-1.))) IF(WDEP.EQ.'y'.OR.WDEP.EQ.'Y')THEN VS=VS0+0.005*(PSI(K,2)+PSI(K-1,2))*DPH1I VSB=VS0+0.01*PSI(K,2)*DPH1I US=US0+U1(K,1) USB=0.5*(U1(K,1)+U1(K+1,1)) VSURF=SQRT(WS0+WD*WD+VS*VS+US*US) VSURFB=SQRT(WS0+WD*WD+VSB*VSB+USB*USB) VSPRIME=VSURF-SQRT(WS0+VS*VS+US*US) ELSE VSURF=SQRT(WS0+WD*WD+VS0*VS0+US0*US0) VSPRIME=VSURF-SQRT(WS0+VS0*VS0+US0*US0) VSURFB=VSURF END IF C FTSURF=G*ROWS*CD*(VSURF*(TS1(K)-TSA)-VSPRIME*TPRIME) 1 *DPH1I FTNET(K,1)=FTNET(K,1)+FTSURF FRSURF=G*ROWS*CD*(VSURF*(RSS-R1(K,1))-VSPRIME*QPRIME) 1 *EFRAC(K)*DPH1I FRNET(K,1)=FRNET(K,1)+FRSURF C FSOIL(K)=0.0 IF(ISOIL.EQ.'y')THEN FRSURFM=100.*3600.*24.*FRSURF/(G*DPH1I) FSOIL(K)=TAUSOILI*(PRECIP+PRADJ-FRSURFM)-TAURUNOFFI*SM1(K) END IF FS=100.*ROWS*CD*CPD*(VSURF*(TS1(K)-TSA)-VSPRIME*TPRIME) FL=100.*2.5E6*ROWS*CD*(VSURF*(RSS-R1(K,1))-VSPRIME*QPRIME) FL=FL*EFRAC(K) C C Put values into time arrays C IF(PRINTIME.GT.(PRINFREQ-10.0))THEN ANPRI=1./FLOAT(NPRINT) ANPRIRADI=1./FLOAT(NPRIRAD) PGRAPH(K,NPR)=PREBAR(K)*ANPRI CAPEGRAPH(K,NPR)=287.*CAPEBAR(K)*ANPRI QRGRAPH(K,NPR)=1000.*QRBAR(K)*ANPRI OGRAPH(K,NPR)=OBAR(K)*ANPRI VGRAPH(K,NPR)=VBAR(K)*ANPRI UTGRAPH(K,NPR)=UBAR(K)*ANPRI SSTGRAPH(K,NPR)=SSTBAR(K)*ANPRI EFRACGRAPH(K,NPR)=EFBAR(K)*ANPRI T100GRAPH(K,NPR)=T100BAR(K)*ANPRI SWGRAPH(K,NPR)=SWBAR(K)*ANPRIRADI OLRGRAPH(K,NPR)=OLRBAR(K)*ANPRIRADI HCGRAPH(K,NPR)=HCBAR(K)*ANPRIRADI LCGRAPH(K,NPR)=LCBAR(K)*ANPRIRADI HCGRAPH(K,NPR)=MIN(HCGRAPH(K,NPR),1.0) LCGRAPH(K,NPR)=MIN(LCGRAPH(K,NPR),1.0) C PREBAR(K)=0.0 CAPEBAR(K)=0.0 QRBAR(K)=0.0 OBAR(K)=0.0 VBAR(K)=0.0 UBAR(K)=0.0 SSTBAR(K)=0.0 EFBAR(K)=0.0 T100BAR(K)=0.0 SWBAR(K)=0.0 OLRBAR(K)=0.0 HCBAR(K)=0.0 LCBAR(K)=0.0 C SST=TS2(K)-273.15 VMAXK(K)=0.0 MPIGRAPH(K,NPR)=0.0 IF(K.LT.KSOIL)THEN CALL PCMIN(SST,PG,P,TCONV,RCONV,NA,N,PMIN,VMAX,IFL) MPIGRAPH(K,NPR)=VMAX VMAXK(K)=VMAX END IF C IF(K.EQ.(NLO-1))THEN GTIME(NPR)=TIME/(3600.*24.) NPR=NPR+1 PRINTIME=0.0 NPRINT=0 NPRIRAD=0 END IF END IF C C Calculate forcing of surface temperature C IF(TSINT.EQ.'y'.OR.TSINT.EQ.'Y')THEN TC=TS1(K)-273.15 ALV=(LV0-CPVMCL*TC)*EFRAC(K) SFLUX=(CPD*(1.-R1(K,1))+CPV*R1(K,1))*100.0*ROWS*CD*(VSURF* 1 (TS1(K)-TSA)-VSPRIME*TPRIME)+ 2 ALV*100.0*ROWS*CD*(VSURF*(RSS-R1(K,1))-VSPRIME*QPRIME) C C Calculate entrainment and ocean current speed and add to ocean heat budget C DELSST=(TS1(K)-TS1(NLO)) DELSST=MAX(DELSST,0.0) C WEH=CCOUP*1.0E-7*SQRT(10.0/(0.1+DELSST))*(VMAXK(K)/70.)**6 C IF(K.LT.(NLO-1))THEN VOCEAN(K)=(VOCEAN(K-1)*COSLAT(K-1)+WEH*DX*COSLATBAR)/COSLAT(K) END IF C FTS(K)=FTS(K)-DELSST*WEH+VOCEAN(K-1)*(TS1(K-1)-TS1(K))/DX C OFLUX(K)=ROWL*4190.0*DEPTH(K)*VOCEAN(K)*DELSST*COSLATBAR C C If radiative cooling is specified, calculate net surface radiative C flux that is consistent with vertically integrated cooling C IF(RADINT.EQ.'y')THEN FTS(K)=FTS(K)+(RADFLUX(K)-SFLUX)/(DEPTH(K)*ROWL*CL) ELSE RADFLC=0.0 DO I=1,NP RADFLC=RADFLC+(CPD*(1.-RCONV(I))+CPV*RCONV(I))*RCOOL(K,I)* 1 100.*(PH(I)-PH(I+1))/G END DO FTS(K)=FTS(K)+(RADFLC-SFLUX)/(DEPTH(K)*ROWL*CL) END IF END IF C C Calculate atmospheric forcings and add to net tendencies C C Level 1 forcings C C Advections: C VOMEG=(COSLAT(K)*PSI(K,2)-COSLAT(K-1)*PSI(K-1,2))*DXI*COSLATBARI+ 1 AOMEGA(2) VVP=0.01*MAX(VOMEG,0.0) VAM=0.01*PSI(K-1,2)*DPH1I+VZ(1) VAP=0.01*PSI(K,2)*DPH1I+VZ(1) VAM=MAX(VAM,0.0) VAP=MIN(VAP,0.0) VA=0.005*(PSI(K-1,2)+PSI(K,2))*DPH1I ALPHA=RD*T2(K,1)*(1.+R2(K,1)*(EPSI-1.))/P(1) C FTADV=VVP*ALPHA/CPD-VVP*(T2(K,1)-T2(K,2))*DP1I FRADV=-VVP*(R2(K,1)-R2(K,2))*DP1I FUADV=-VVP*(U2(K,1)-U2(K,2))*DP1I FTHADV=-(VAM*(T2(K,1)-T2(K-1,1))+VAP*(T2(K+1,1)-T2(K,1)))*DXI FRHADV=-(VAM*(R2(K,1)-R2(K-1,1))+VAP*(R2(K+1,1)-R2(K,1)))*DXI FUHADV=-(VAM*(U2(K,1)-U2(K-1,1))+VAP*(U2(K+1,1)-U2(K,1)))*DXI C AFLUX(K)=PSI(K,2)*0.5*(CPD*(T2(K+1,1)+T2(K,1))+LV0* C 1 (R2(K+1,1)+R2(K,1))) AFLUX(K)=-(CPD*(FTHADV+FTADV)+LV0*(FRHADV+FRADV))* 1 (PH(1)-PH(2)) C UF(K,1)=VA*FCBAR+FUADV+FUHADV+FUC(1) c UF(K,1)=0.0 ETAF(K,1)=0.0 FTUADV=-U2(K,1)*TX(1) FRUADV=5423.0*R1(K,1)*FTUADV/(T2(K,1)*T2(K,1)) FTNET(K,1)=FTNET(K,1)+FTADV+FTHADV+FTUADV FRNET(K,1)=FRNET(K,1)+FRADV+FRHADV+FRUADV C C Level 1 Vertical diffusion terms C PNUP=PNUM*(P(2)-P(1)+DPBL)/DPBL PNUP=MAX(PNUP,0.0) c TAUETA0=CD*VSURFB*ETA1(K,1)*ROWS*G TAUETA0=-CD*VSURFB*PSI(K,2)*1.0E-4*ROWS*G*DPH1I*DPH1I DIFETA=2.*PNUP*(ETA1(K,2)-ETA1(K,1))-TAUETA0*DPH1I C PNUP=PNUM*(PH(2)-PH(1)+DPBL)/DPBL PNUP=MAX(PNUP,0.0) TAUU0=CD*VSURF*(U1(K,1)+US0)*ROWS*G DIFU=2.*PNUP*(U1(K,2)-U1(K,1))-TAUU0*DP1I C ETAF(K,1)=ETAF(K,1)+DIFETA ETAF(K,2)=DIFETA UF(K,1)=UF(K,1)+DIFU C C Level 1 Horizontal diffusion terms C DIFHU=2.*DAMP*(U1(K+1,1)*COSLAT(K+1)+U1(K-1,1)*COSLAT(K-1)- 1 2.*U1(K,1)*COSLAT(K))/COSLAT(K) DIFHT=2.*DAMP*(T1(K+1,1)+T1(K-1,1)-2.*T1(K,1)) DIFHR=2.*DAMP*(R1(K+1,1)+R1(K-1,1)-2.*R1(K,1)) c DIFHR=0.0 c DIFHT=0.0 AFLUX(K)=AFLUX(K)-DIFHT*CPD*(PH(1)-PH(2)) C UF(K,1)=UF(K,1)+DIFHU FTNET(K,1)=FTNET(K,1)+DIFHT FRNET(K,1)=FRNET(K,1)+DIFHR GZ1=0.0 GZ2=0.0 C C Interior forcings C DO 140 I=2, NP-1 C C Advections: C VOMEG1=(COSLAT(K)*PSI(K,I+1)-COSLAT(K-1)*PSI(K-1,I+1))*DXI* 1 COSLATBARI+AOMEGA(I+1) VOMEG=(COSLAT(K)*PSI(K,I)-COSLAT(K-1)*PSI(K-1,I))*DXI* 1 COSLATBARI+AOMEGA(I) VVP=0.01*MAX(VOMEG1,0.0) VVM=0.01*MIN(VOMEG,0.0) ALPHA=RD*T2(K,I)*(1.+R2(K,I)*(EPSI-1.))/P(I) VAM=0.01*(PSI(K-1,I+1)-PSI(K-1,I))/(PH(I)-PH(I+1))+VZ(I) VAP=0.01*(PSI(K,I+1)-PSI(K,I))/(PH(I)-PH(I+1))+VZ(I) VA=0.5*(VAM+VAP)-VZ(I) VAM=MAX(VAM,0.0) VAP=MIN(VAP,0.0) C VA=VAP+VAM FTADV=(VVP+VVM)*ALPHA/CPD- 1 VVP*(T2(K,I)-T2(K,I+1))/(P(I)-P(I+1))- 2 VVM*(T2(K,I-1)-T2(K,I))/(P(I-1)-P(I)) FRADV=-VVP*(R2(K,I)-R2(K,I+1))/(P(I)-P(I+1))- 1 VVM*(R2(K,I-1)-R2(K,I))/(P(I-1)-P(I)) FUADV=-VVP*(U2(K,I)-U2(K,I+1))/(P(I)-P(I+1))- 1 VVM*(U2(K,I-1)-U2(K,I))/(P(I-1)-P(I)) FTHADV=-(VAM*(T2(K,I)-T2(K-1,I))+VAP*(T2(K+1,I)-T2(K,I)))*DXI FRHADV=-(VAM*(R2(K,I)-R2(K-1,I))+VAP*(R2(K+1,I)-R2(K,I)))*DXI FUHADV=-(VAM*(U2(K,I)-U2(K-1,I))+VAP*(U2(K+1,I)-U2(K,I)))*DXI C UF(K,I)=VA*FCBAR+FUADV+FUHADV+FUC(I) c UF(K,I)=0.0 FTUADV=-U2(K,I)*TX(I) FRUADV=5423.0*R1(K,I)*FTUADV/(T2(K,I)*T2(K,I)) C FTNET(K,I)=FTNET(K,I)+FTADV+FTHADV+FTUADV FRNET(K,I)=FRNET(K,I)+FRADV+FRHADV+FRUADV AFLUX(K)=AFLUX(K)-(CPD*(FTHADV+FTADV)+LV0*(FRHADV+FRADV))* 1 (PH(I)-PH(I+1)) C C Calculate forcing of horizontal vorticity C TVDIF1=T2(K+1,I)*(1.+R2(K+1,I)*(EPSI-1.))- 1 T2(K,I)*(1.+R2(K,I)*(EPSI-1.)) TVDIF2=T2(K+1,I-1)*(1.+R2(K+1,I-1)*(EPSI-1.))- 1 T2(K,I-1)*(1.+R2(K,I-1)*(EPSI-1.)) ETAF(K,I)=0.005*RD*DXI*(TVDIF1/PH(I)+TVDIF2/PH(I-1)) ETAF(K,I)=ETAF(K,I)-FC(K)*0.005*(U2(K,I-1)+U2(K+1,I-1)- 1 U2(K,I)-U2(K+1,I))/(P(I-1)-P(I)) C C Convective momentum flux forcing of vorticity C ETAF(K-1,I)=ETAF(K-1,I)+0.005*(FVCNET(K,I-1)-FVCNET(K,I)+ 1 FVCNET(K-1,I-1)-FVCNET(K-1,I))/(PH(I-1)-PH(I)) C C Advection of vorticity by horizontal flow C VAM=0.5*(VZ(I-1)+VZ(I))+0.005*(PSI(K,I+1)-PSI(K,I-1)+ 1 PSI(K-1,I+1)-PSI(K-1,I-1))/(PH(I-1)-PH(I+1)) VAP=0.5*(VZ(I-1)+VZ(I))+0.005*(PSI(K,I+1)-PSI(K,I-1)+ 1 PSI(K+1,I+1)-PSI(K+1,I-1))/(PH(I-1)-PH(I+1)) VAM=MAX(VAM,0.0) VAP=MIN(VAP,0.0) ETAF(K,I)=ETAF(K,I)-(VAM*(ETA2(K,I)- 1 ETA2(K-1,I))+VAP*(ETA2(K+1,I)-ETA2(K,I)))*DXI C C Vertical diffusion terms C PNU=PNUM*(P(I)-P(1)+DPBL)/DPBL PNU=MAX(PNU,0.0) PNUP=PNUM*(P(I+1)-P(1)+DPBL)/DPBL PNUP=MAX(PNUP,0.0) DIFETA=(PH(I-1)-PH(I+1))*(PNUP*(ETA1(K,I+1)-ETA1(K,I))/ 1 (PH(I)-PH(I+1))-PNU*(ETA1(K,I)-ETA1(K,I-1))/ 2 (PH(I-1)-PH(I))) C PNU=PNUM*(PH(I)-PH(1)+DPBL)/DPBL PNU=MAX(PNU,0.0) PNUP=PNUM*(PH(I+1)-PH(1)+DPBL)/DPBL PNUP=MAX(PNUP,0.0) DIFU=(P(I-1)-P(I+1))*(PNUP*(U1(K,I+1)-U1(K,I))/ 1 (P(I)-P(I+1))-PNU*(U1(K,I)-U1(K,I-1))/(P(I-1)-P(I))) c DIFU=10.*DIFU C ETAF(K,I)=ETAF(K,I)+DIFETA UF(K,I)=UF(K,I)+DIFU C C Horizontal diffusion terms C DIFHETA=2.*DAMP*(ETA1(K+1,I)+ETA1(K-1,I)-2.*ETA1(K,I)) DIFHU=2.*DAMP*(U1(K+1,I)*COSLAT(K+1)+U1(K-1,I)*COSLAT(K-1)- 1 2.*U1(K,I)*COSLAT(K))/COSLAT(K) DIFHT=2.*DAMP*(T1(K+1,I)+T1(K-1,I)-2.*T1(K,I)) DIFHR=2.*DAMP*(R1(K+1,I)+R1(K-1,I)-2.*R1(K,I)) c DIFHT=0.0 c DIFHR=0.0 AFLUX(K)=AFLUX(K)-DIFHT*CPD*(PH(I)-PH(I+1)) C ETAF(K,I)=ETAF(K,I)+DIFHETA UF(K,I)=UF(K,I)+DIFHU FTNET(K,I)=FTNET(K,I)+DIFHT FRNET(K,I)=FRNET(K,I)+DIFHR C C Sponge layer damping C FTNET(K,I)=FTNET(K,I)+SPONGE(I)*(T0(K,I)-T1(K,I)) FRNET(K,I)=FRNET(K,I)+SPONGE(I)*(R0(K,I)-R1(K,I)) UF(K,I)=UF(K,I)-SPONGE(I)*U1(K,I) ETAF(K,I)=ETAF(K,I)-SPONGE(I)*ETA1(K,I) C C IF(I.GE.NPTOP)THEN C ETAF(K,I)=0.0 C UF(K,I)=0.0 C END IF C 140 CONTINUE AFLUX(K)=AFLUX(K-1)+AFLUX(K)*100.*DX*COSLATBAR/G C C Level NP forcings C UF(K,NP)=0.0 C C Calculate certain graphics quantities C DO 150 I=1,NP MUPG(K,I)=1000.*MUP(I) MDG(K,I)=-1000.*MP(I) RHG(K,I)=100.*R2(K,I)/(RS(I)+2.0E-6) RHG(K,I)=MIN(RHG(K,I),100.0) RHG(K,I)=MAX(RHG(K,I),0.0) OLG(K,I)=36.*(COSLAT(K)*PSI(K,I)-COSLAT(K-1)*PSI(K-1,I))* 1 DXI*COSLATBARI 150 CONTINUE C 155 CONTINUE OFLUX(NLO-1)=0.0 AFLUX(NLO-1)=0.0 C DO I=2,NP-1 ETAF(NLO-1,I)=ETAF(NLO-1,I)+0.005*(FVCNET(NLO-1,I-1)- 1 FVCNET(NLO-1,I)+FVCNET(NLO-2,I-1)-FVCNET(NLO-2,I))/ 2 (PH(I-1)-PH(I)) END DO C C Calculate moving average plots C IF(TIME.GT.(ENDTIME-PAVE))THEN AFLEN=2.*3.14159*6.38E6*1.0E-15*PAFRAC DO K=2,NLO-1 AFLUXGRAPH(K)=AFLUXGRAPH(K)+AFLUX(K)*AFLEN OFLUXGRAPH(K)=OFLUXGRAPH(K)+OFLUX(K)*AFLEN NETFLUXGRAPH(K)=NETFLUXGRAPH(K)+NETFLUX(K)*AFLEN END DO DO 203 I=1,NP DO 203 K=2,NLO-1 J=K+(VTRANS*TIME*DXI) IF(VTRANS.GE.0.0)THEN J=J-(NLO-2)*(J/NLO) ELSE J=J+(NLO-1)*((NLO-J)/NLO) END IF OLGRAPH(K,I)=OLGRAPH(K,I)+OLG(J,I)*PAFRAC MUPGRAPH(K,I)=MUPGRAPH(K,I)+MUPG(J,I)*PAFRAC MDGRAPH(K,I)=MDGRAPH(K,I)+MDG(J,I)*PAFRAC RHGRAPH(K,I)=RHGRAPH(K,I)+RHG(J,I)*PAFRAC ULGRAPH(K,I)=ULGRAPH(K,I)+U2(J,I)*PAFRAC PSIGRAPH(K,I)=PSIGRAPH(K,I)+1.0E-5*PSI(J,I)* 1 COS(PIM*RLAT(J))*PAFRAC TEMPGRAPH(K,I)=TEMPGRAPH(K,I)+T2(J,I)*PAFRAC TEMPGRAPHNET(K,I)=TEMPGRAPH(K,I) AMGRAPH(K,I)=AMGRAPH(K,I)+(U2(J,I)+ 1 1.454E-4*6.38E6*COS(PIM*0.5*(RLAT(J)+RLAT(J-1))))* 2 PAFRAC CLDGRAPH(K,I)=CLDGRAPH(K,I)+CLDTEMP(K,I)*PAFRAC 203 CONTINUE END IF C C Advance quantities one time step, applying Asselin filter C DO 180 K=2,NLO-1 C IF(TSINT.EQ.'y'.OR.TSINT.EQ.'Y')THEN TS3(K)=TS1(K)+2.*DT*FTS(K) TS3(K)=MAX(TS3(K),270.0) TS1(K)=TS2(K)+AFILT*(TS1(K)+TS3(K)-2.*TS2(K)) TS2(K)=TS3(K) END IF C SM3(K)=SM1(K)+2.*DT*FSOIL(K) SM1(K)=SM2(K)+AFILT*(SM1(K)+SM3(K)-2.*SM2(K)) SM2(K)=SM3(K) C DO 160 I=1,NP T3(K,I)=T1(K,I)+2.*DT*FTNET(K,I) R3(K,I)=R1(K,I)+2.*DT*FRNET(K,I) T1(K,I)=T2(K,I)+AFILT*(T1(K,I)+T3(K,I)-2.*T2(K,I)) R1(K,I)=R2(K,I)+AFILT*(R1(K,I)+R3(K,I)-2.*R2(K,I)) T2(K,I)=T3(K,I) R2(K,I)=R3(K,I) IF(I.LE.NPTOP)THEN ETA3(K,I)=ETA1(K,I)+2.*DT*ETAF(K,I) ETA1(K,I)=ETA2(K,I)+AFILT*(ETA1(K,I)+ETA3(K,I)-2.*ETA2(K,I)) ETA2(K,I)=ETA3(K,I) U3(K,I)=U1(K,I)+2.*DT*UF(K,I) U1(K,I)=U2(K,I)+AFILT*(U1(K,I)+U3(K,I)-2.*U2(K,I)) U2(K,I)=U3(K,I) END IF 160 CONTINUE ETA1(K,1)=ETA1(K,2) C C Find streamfunction C PSI(K,NPTOP)=0.0 ETA2(K,NPTOP)=0.0 FV=0.0 DO 165 I=NPTOP-1,1,-1 FV=FV+ETA2(K,I+1)*100.0*(P(I)-P(I+1)) PSI(K,I)=PSI(K,I+1)+FV*100.0*(PH(I)-PH(I+1)) 165 CONTINUE DO 170 I=NPTOP-1,2,-1 PSI(K,I)=PSI(K,1)*(PH(I)-PH(NPTOP))/(PH(1)-PH(NPTOP))-PSI(K,I) 170 CONTINUE PSI(K,1)=0.0 C 180 CONTINUE C C Return to beginning of time loop C GOTO 100 C C Output processing starts here C 200 CONTINUE NPR=NPR-1 C IF(FWRITE.EQ.'y')THEN FNI=1./FLOAT(NSAVE) DO I=1,NP RV2M(I)=FNI*RV2M(I) CLDFRADM(I)=FNI*CLDFRADM(I) CLDEMIM(I)=FNI*CLDEMIM(I) CLDTAUM(I)=FNI*CLDTAUM(I) END DO END IF C C Subtract horizontal mean temperature from composite temperature C DO 1010 I=1,NP II=MIN(I,NP-1) TEMBAR=0.0 DO 1005 K=2,NLO-1 TEMBAR=TEMBAR+TEMPGRAPH(K,II) 1005 CONTINUE TEMBAR=TEMBAR/FLOAT(NLO-2) DO 1006 K=2,NLO-1 TEMPGRAPH(K,I)=TEMPGRAPH(K,II)-TEMBAR 1006 CONTINUE 1010 CONTINUE C C *** Add boundary conditions to contour plots *** C IF(BCS.EQ.'p')THEN C DO 205 I=1,NPR PGRAPH(1,I)=PGRAPH(NLO-1,I) PGRAPH(NLO,I)=PGRAPH(2,I) CAPEGRAPH(1,I)=CAPEGRAPH(NLO-1,I) CAPEGRAPH(NLO,I)=CAPEGRAPH(2,I) QRGRAPH(1,I)=QRGRAPH(NLO-1,I) QRGRAPH(NLO,I)=QRGRAPH(2,I) OGRAPH(1,I)=OGRAPH(NLO-1,I) OGRAPH(NLO,I)=OGRAPH(2,I) VGRAPH(1,I)=VGRAPH(NLO-1,I) VGRAPH(NLO,I)=VGRAPH(2,I) UTGRAPH(1,I)=UTGRAPH(NLO-1,I) UTGRAPH(NLO,I)=UTGRAPH(2,I) SSTGRAPH(1,I)=SSTGRAPH(NLO-1,I) SSTGRAPH(NLO,I)=SSTGRAPH(2,I) EFRACGRAPH(1,I)=EFRACGRAPH(NLO-1,I) EFRACGRAPH(NLO,I)=EFRACGRAPH(2,I) T100GRAPH(1,I)=T100GRAPH(NLO-1,I) T100GRAPH(NLO,I)=T100GRAPH(2,I) MPIGRAPH(1,I)=MPIGRAPH(NLO-1,I) MPIGRAPH(NLO,I)=MPIGRAPH(2,I) SWGRAPH(1,I)=SWGRAPH(NLO-1,I) SWGRAPH(NLO,I)=SWGRAPH(2,I) OLRGRAPH(1,I)=OLRGRAPH(NLO-1,I) OLRGRAPH(NLO,I)=OLRGRAPH(2,I) HCGRAPH(1,I)=HCGRAPH(NLO-1,I) HCGRAPH(NLO,I)=HCGRAPH(2,I) LCGRAPH(1,I)=LCGRAPH(NLO-1,I) LCGRAPH(NLO,I)=LCGRAPH(2,I) 205 CONTINUE DO 207 J=1,NP MUPGRAPH(1,J)=MUPGRAPH(NLO-1,J) MUPGRAPH(NLO,J)=MUPGRAPH(2,J) MDGRAPH(1,J)=MDGRAPH(NLO-1,J) MDGRAPH(NLO,J)=MDGRAPH(2,J) OLGRAPH(1,J)=OLGRAPH(NLO-1,J) OLGRAPH(NLO,J)=OLGRAPH(2,J) RHGRAPH(1,J)=RHGRAPH(NLO-1,J) RHGRAPH(NLO,J)=RHGRAPH(2,J) ULGRAPH(1,J)=ULGRAPH(NLO-1,J) ULGRAPH(NLO,J)=ULGRAPH(2,J) PSIGRAPH(1,J)=PSIGRAPH(NLO-1,J) PSIGRAPH(NLO,J)=PSIGRAPH(2,J) TEMPGRAPH(1,J)=TEMPGRAPH(NLO-1,J) TEMPGRAPH(NLO,J)=TEMPGRAPH(2,J) TEMPGRAPHNET(1,J)=TEMPGRAPHNET(NLO-1,J) TEMPGRAPHNET(NLO,J)=TEMPGRAPHNET(2,J) AMGRAPH(1,J)=AMGRAPH(NLO-1,J) AMGRAPH(NLO,J)=AMGRAPH(2,J) CLDGRAPH(1,J)=CLDGRAPH(NLO-1,J) CLDGRAPH(NLO,J)=CLDGRAPH(2,J) 207 CONTINUE C ELSE C DO 206 I=1,NPR PGRAPH(1,I)=PGRAPH(2,I) PGRAPH(NLO,I)=PGRAPH(NLO-1,I) CAPEGRAPH(1,I)=CAPEGRAPH(2,I) CAPEGRAPH(NLO,I)=CAPEGRAPH(NLO-1,I) QRGRAPH(1,I)=QRGRAPH(2,I) QRGRAPH(NLO,I)=QRGRAPH(NLO-1,I) OGRAPH(1,I)=OGRAPH(2,I) OGRAPH(NLO,I)=OGRAPH(NLO-1,I) VGRAPH(1,I)=0.0 VGRAPH(NLO,I)=0.0 UTGRAPH(1,I)=0.0 UTGRAPH(NLO,I)=0.0 SSTGRAPH(1,I)=SSTGRAPH(2,I) SSTGRAPH(NLO,I)=SSTGRAPH(NLO-1,I) EFRACGRAPH(1,I)=EFRACGRAPH(2,I) EFRACGRAPH(NLO,I)=EFRACGRAPH(NLO-1,I) T100GRAPH(1,I)=T100GRAPH(2,I) T100GRAPH(NLO,I)=T100GRAPH(NLO-1,I) MPIGRAPH(1,I)=MPIGRAPH(2,I) MPIGRAPH(NLO,I)=MPIGRAPH(NLO-1,I) SWGRAPH(1,I)=SWGRAPH(2,I) SWGRAPH(NLO,I)=SWGRAPH(NLO-1,I) OLRGRAPH(1,I)=OLRGRAPH(2,I) OLRGRAPH(NLO,I)=OLRGRAPH(NLO-1,I) HCGRAPH(1,I)=HCGRAPH(2,I) HCGRAPH(NLO,I)=HCGRAPH(NLO-1,I) LCGRAPH(1,I)=LCGRAPH(2,I) LCGRAPH(NLO,I)=LCGRAPH(NLO-1,I) 206 CONTINUE DO 208 J=1,NP MUPGRAPH(1,J)=MUPGRAPH(2,J) MUPGRAPH(NLO,J)=MUPGRAPH(NLO-1,J) MDGRAPH(1,J)=MDGRAPH(2,J) MDGRAPH(NLO,J)=MDGRAPH(NLO-1,J) OLGRAPH(1,J)=OLGRAPH(2,J) OLGRAPH(NLO,J)=OLGRAPH(NLO-1,J) RHGRAPH(1,J)=RHGRAPH(2,J) RHGRAPH(NLO,J)=RHGRAPH(NLO-1,J) ULGRAPH(1,J)=0.0 ULGRAPH(NLO,J)=0.0 PSIGRAPH(1,J)=0.0 PSIGRAPH(NLO,J)=0.0 TEMPGRAPH(1,J)=TEMPGRAPH(2,J) TEMPGRAPH(NLO,J)=TEMPGRAPH(NLO-1,J) TEMPGRAPHNET(1,J)=TEMPGRAPHNET(2,J) TEMPGRAPHNET(NLO,J)=TEMPGRAPHNET(NLO-1,J) AMGRAPH(NLO,J)=AMGRAPH(NLO-1,J) AMGRAPH(1,J)=AMGRAPH(2,J) CLDGRAPH(1,J)=CLDGRAPH(2,J) CLDGRAPH(NLO,J)=CLDGRAPH(NLO-1,J) 208 CONTINUE C END IF C C *** Write contour plot ASCII files *** C DO 400 I=1,NLO DISX(I)=0.001*DX*FLOAT(I-1) 400 CONTINUE C OPEN(UNIT=12, FILE='output/x.out',STATUS='UNKNOWN') WRITE(12,425)(DISX(I),I=1,NLO) CLOSE(12) C OPEN(UNIT=12, FILE='output/x2.out',STATUS='UNKNOWN') WRITE(12,425)((DISX(I)+0.0005*DX),I=1,NLO-1) CLOSE(12) C OPEN(UNIT=12, FILE='output/t.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,435)GTIME(I) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/p.out',STATUS='UNKNOWN') WRITE(12,437)(P(I),I=1,NP) CLOSE(12) C 425 FORMAT(2X,200(F9.1,' ')) 435 FORMAT(2X,500(F14.5,' ')) 437 FORMAT(2X,100(F9.1,' ')) C OPEN(UNIT=12, FILE='output/precip.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(PGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/omegat.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(OGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/vt.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(VGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/ut.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(UTGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/sst.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(SSTGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/efrac.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(EFRACGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/t100.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(T100GRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/mpi.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(MPIGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/sw.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(SWGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/olr.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(OLRGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/hc.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(HCGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/lc.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(LCGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12,FILE='output/cape.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(CAPEGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12,FILE='output/qr.out',STATUS='UNKNOWN') DO I=1,NPR WRITE(12,445)(QRGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/omegal.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(OLGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/mul.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(MUPGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/mdl.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(MDGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/psi.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(PSIGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/rh.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(RHGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/cld.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(CLDGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/ul.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(ULGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/temp.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(TEMPGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/am.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(AMGRAPH(J,I),J=1,NLO) END DO CLOSE(12) C OPEN(UNIT=12, FILE='output/tempnet.out',STATUS='UNKNOWN') DO I=1,NP WRITE(12,445)(TEMPGRAPHNET(J,I),J=1,NLO) END DO CLOSE(12) C 445 FORMAT(2X,200(F10.4,' ')) 450 CONTINUE C OPEN(UNIT=12, FILE='output/flux.out',STATUS='UNKNOWN') DO K=1,NLO-1 WRITE(12,*)(DISX(K)+0.0005*DX),AFLUXGRAPH(K),OFLUXGRAPH(K), 1 NETFLUXGRAPH(K) END DO CLOSE(12) C C Output to file sounding.out C OPEN(UNIT=10, FILE='output/sounding.out', STATUS='unknown') C DO 350 K=2,NLO-1 TS2C=TS2(K)-273.15 WRITE(10,210) NP,QGRAPH(K,NPR),TS2C 210 FORMAT(2X, 'N levels = ',T15,I3,T26,'Q=',F5.2,' m/s', 1 T48,'SST=',T55,F5.2,' C',/) C WRITE(10,220) 220 FORMAT(2X,'Pressure',T14,'Temp.',T25,'Spec. humid.',T40, 1 'O3 mr',T49, 'Psi',T57,'Rcool',T67,'U') WRITE(10,230) 230 FORMAT(4X,'(mb)',T15,'(C)',T28,'(g/kg)',T40,'10**-6',T47, 1 'E4 kg/m/s',T57,'(K/dy)',T66,'m/s') WRITE(10,240) 240 FORMAT(2X,'-------',T14,'-----',T25,'------------',T40, 1 '-----',T48, '-----',T57,'------',T67,'-') c DO 260 I=1,NP OZ=1.0E6*O3(I) TC=T2(K,I)-273.15 RM=1000.0*R2(K,I) AOMEG=PSI(K,I)*1.0E-4 ARADC=RCOOL(K,I)*3600.*24. WRITE(10,250)P(I),TC,RM,OZ,AOMEG,ARADC,U2(K,I) 250 FORMAT(3X,F6.1,T12,F8.3,T25,F10.6,T39,F6.3,T47,F6.2, 1 T56,F7.3,T64,F5.1) 260 CONTINUE 350 CONTINUE CLOSE(10) IF(FWRITE.EQ.'y')THEN OPEN(10,FILE='output/variables',form='unformatted', 1 status='unknown') WRITE(10) RV2M,CLDFRADM,CLDEMIM,CLDTAUM CLOSE(10) END IF STOP END C SUBROUTINE PCMIN(SST,PSL,P,T,R,NA,N,PMIN,VMAX,IFL) C C *** This subroutine calculates the mimimum central pressure *** C *** and maximum 10 m 1 minute wind speed *** C *** achievable in tropical cyclones, given a sounding *** C *** and a sea surface temperature. *** C C INPUT: SST: Sea surface temperature in C C C PSL: Sea level pressure (mb) C C P,T,R: One-dimensional arrays of dimension NA C containing pressure (mb), temperature (K), C and mixing ratio (kg/kg). The arrays MUST be C arranged so that the lowest index corresponds C to the lowest model level, with increasing index C corresponding to decreasing pressure. The temperature C sounding should extend to at least the tropopause and C preferably to the lower stratosphere, however the C mixing ratios are not important above the boundary C layer. Missing mixing ratios can be replaced by zeros. C C NA: The dimension of P,T and R C C N: The actual number of points in the sounding C (N is less than or equal to NA) C C OUTPUT: PMIN is the minimum central pressure, in mb C C VMAX is the maximum surface wind speed, in m/s C (reduced to reflect surface drag) C C IFL is a flag: A value of 1 means OK; a value of 0 C indicates no convergence (hypercane); a value of 2 C means that the CAPE routine failed. C C----------------------------------------------------------------------------- REAL T(NA), P(NA), R(NA) C C *** Adjustable constant: Ratio of C_k to C_D *** C CKCD=0.9 C C *** Adjustable constant for buoyancy of displaced parcels: *** C *** 0=Reversible ascent; 1=Pseudo-adiabatic ascent *** C SIG=0.5 C C *** Adjustable switch: if IDISS = 0, no dissipative heating is *** C *** allowed; otherwise, it is *** C IDISS=1 C C *** Normalize certain quantities *** C SSTK=SST+273.15 ES0=6.112*EXP(17.67*SST/(243.5+SST)) C C *** Default values *** C VMAX=0.0 PMIN=PSL IFL=1 C NP=0 PM=950.0 C C *** Find environmental CAPE *** C TP=T(1) RP=R(1) PP=P(1) CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEA,TOA,IFLAG) IF(IFLAG.NE.1)IFL=2 C C *** Begin iteration to find mimimum pressure *** C 100 CONTINUE C C *** Find CAPE at radius of maximum winds *** C TP=T(1) PP=PM RP=0.622*R(1)*PSL/(PM*(0.622+R(1))-R(1)*PSL) CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEM,TOM,IFLAG) IF(IFLAG.NE.1)IFL=2 RAT=SSTK/MAX(TOM,1.0) IF(IDISS.EQ.0)RAT=1.0 C C *** Find saturation CAPE at radius of maximum winds *** C TP=SSTK PP=PM RP=0.622*ES0/(PM-ES0) CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEMS,TOMS,IFLAG) IF(IFLAG.NE.1)IFL=2 C C *** Initial estimate of minimum pressure *** C RS0=RP TV1=T(1)*(1.+R(1)/0.622)/(1.+R(1)) TVAV=0.5*(TV1+SSTK*(1.+RS0/0.622)/(1.+RS0)) CAT=CAPEM-CAPEA+0.5*CKCD*RAT*(CAPEMS-CAPEM) PNEW=PSL*EXP(-CAT/(287.04*TVAV)) C C *** Test for convergence *** C IF(ABS(PNEW-PM).GT.0.2)THEN PM=PNEW NP=NP+1 IF(NP.GT.1000.OR.PM.LT.400.0)THEN PMIN=400.0 IFL=0 GOTO 900 END IF GOTO 100 ELSE CAT=CAPEM-CAPEA+CKCD*RAT*(CAPEMS-CAPEM) PMIN=PSL*EXP(-CAT/(287.04*TVAV)) END IF 900 CONTINUE FAC=MAX(0.0,(CAPEMS-CAPEM)) VMAX=0.8*SQRT(CKCD*RAT*FAC) C C *** Renormalize sounding arrays *** C DO 910 I=1,N R(I)=R(I)*1000.0 T(I)=T(I)-273.15 910 CONTINUE C RETURN END C SUBROUTINE CAPE(TP,RP,PP,T,R,P,ND,N,SIG,CAPED,TO,IFLAG) C C This subroutine calculates the CAPE of a parcel with pressure PP (mb), C temperature TP (K) and mixing ratio RP (gm/gm) given a sounding C of temperature (T in K) and mixing ratio (R in gm/gm) as a function C of pressure (P in mb). ND is the dimension of the arrays T,R and P, C while N is the actual number of points in the sounding. CAPED is C the calculated value of CAPE and TO is the temperature at the C level of neutral buoyancy. IFLAG is a flag C integer. If IFLAG = 1, routine is successful; if it is 0, routine did C not run owing to improper sounding (e.g.no water vapor at parcel level). C IFLAG=2 indicates that routine did not converge. C REAL T(ND),R(ND),P(ND),TVRDIF(100) REAL NA C C *** Default values *** C CAPED=0.0 TO=T(1) IFLAG=1 C C *** Check that sounding is suitable *** C IF(RP.LT.1.0E-6.OR.TP.LT.200.0)THEN IFLAG=0 RETURN END IF C C *** Assign values of thermodynamic constants *** C CPD=1005.7 CPV=1870.0 CL=2500.0 C CPVMCL=2320.0 CPVMCL=CL-CPV RV=461.5 RD=287.04 EPS=RD/RV ALV0=2.501E6 C C *** Define various parcel quantities, including reversible *** C *** entropy, S. *** C TPC=TP-273.15 ESP=6.112*EXP(17.67*TPC/(243.5+TPC)) EVP=RP*PP/(EPS+RP) RH=EVP/ESP ALV=ALV0-CPVMCL*TPC S=(CPD+RP*CL)*LOG(TP)-RD*LOG(PP-EVP)+ 1 ALV*RP/TP-RP*RV*LOG(RH) C C *** Find lifted condensation pressure, PLCL *** C CHI=TP/(1669.0-122.0*RH-TP) PLCL=PP*(RH**CHI) C C *** Begin updraft loop *** C NCMAX=0 DO J=1,N TVRDIF(J)=0.0 END DO DO 200 J=2,N C C *** Don't bother lifting parcel above 60 mb *** C IF(P(J).LT.59.0)GOTO 200 C C *** Parcel quantities below lifted condensation level *** C IF(P(J).GE.PLCL)THEN TG=TP*(P(J)/PP)**(RD/CPD) RG=RP C C *** Calculate buoyancy *** C TLVR=TG*(1.+RG/EPS)/(1.+RG) TVRDIF(J)=TLVR-T(J)*(1.+R(J)/EPS)/(1.+R(J)) ELSE C C *** Parcel quantities above lifted condensation level *** C TG=T(J) TJC=T(J)-273.15 ES=6.112*EXP(17.67*TJC/(243.5+TJC)) RG=EPS*ES/(P(J)-ES) C C *** Iteratively calculate lifted parcel temperature and mixing *** C *** ratio for reversible ascent *** C NC=0 120 CONTINUE NC=NC+1 C C *** Calculate estimates of the rates of change of the entropy *** C *** with temperature at constant pressure *** C ALV=ALV0-CPVMCL*(TG-273.15) SL=(CPD+RP*CL+ALV*ALV*RG/(RV*TG*TG))/TG EM=RG*P(J)/(EPS+RG) SG=(CPD+RP*CL)*LOG(TG)-RD*LOG(P(J)-EM)+ 1 ALV*RG/TG IF(NC.LT.3)THEN AP=0.3 ELSE AP=1.0 END IF TGNEW=TG+AP*(S-SG)/SL C C *** Test for convergence *** C IF(ABS(TGNEW-TG).GT.0.01)THEN TG=TGNEW TC=TG-273.15 ENEW=6.112*EXP(17.67*TC/(243.5+TC)) C C *** Bail out if things get out of hand *** C IF(NC.GT.500.OR.ENEW.GT.(P(J)-1.0))THEN IFLAG=2 RETURN END IF RG=EPS*ENEW/(P(J)-ENEW) GOTO 120 END IF NCMAX=MAX(NC,NCMAX) C C *** Calculate buoyancy *** C RMEAN=SIG*RG+(1.-SIG)*RP TLVR=TG*(1.+RG/EPS)/(1.+RMEAN) TVRDIF(J)=TLVR-T(J)*(1.+R(J)/EPS)/(1.+R(J)) END IF 200 CONTINUE C C *** Begin loop to find NA, PA, and CAPE from reversible ascent *** C NA=0.0 PA=0.0 C C *** Find maximum level of positive buoyancy, INB *** C INB=1 DO 550 J=N,1,-1 IF(TVRDIF(J).GT.0.0)INB=MAX(INB,J) 550 CONTINUE C C *** Find positive and negative areas and CAPE *** C IF(INB.GT.1)THEN DO 600 J=2,INB TVM=0.5*(TVRDIF(J)+TVRDIF(J-1)) PMA=0.5*(P(J)+P(J-1)) IF(TVM.LE.0.0)THEN NA=NA-RD*TVM*(P(J-1)-P(J))/PMA ELSE PA=PA+RD*TVM*(P(J-1)-P(J))/PMA END IF 600 CONTINUE C C *** Find residual positive area above INB and TO *** C PAT=0.0 TO=T(INB) IF(INB.LT.N)THEN PINB=(P(INB+1)*TVRDIF(INB)-P(INB)*TVRDIF(INB+1))/ 1 (TVRDIF(INB)-TVRDIF(INB+1)) PAT=RD*TVRDIF(INB)*(P(INB)-PINB)/(P(INB)+PINB) TO=(T(INB)*(PINB-P(INB+1))+T(INB+1)*(P(INB)-PINB))/ 1 (P(INB)-P(INB+1)) END IF C C *** Find CAPE *** C CAPED=PA+PAT-NA CAPED=MAX(CAPED,0.0) END IF C RETURN END C *************************************************************************** ***** SUBROUTINE CONVECT ***** ***** VERSION 4.3b ***** ***** 20 August, 2000 ***** ***** Kerry Emanuel ***** *************************************************************************** C SUBROUTINE CONVECT * (T, Q, QS, U, V, TRA, P, PH, * ND, NL, NTRA, DELT, IFLAG, FT, FQ, FU, * FV, FTRA, PRECIP, WD, TPRIME, QPRIME, QCONDC,CBMF,FS,FL) C C----------------------------------------------------------------------------- C *** On input: *** C C T: Array of absolute temperature (K) of dimension ND, with first C index corresponding to lowest model level. Note that this array C will be altered by the subroutine if dry convective adjustment C occurs and if IPBL is not equal to 0. C C Q: Array of specific humidity (gm/gm) of dimension ND, with first C index corresponding to lowest model level. Must be defined C at same grid levels as T. Note that this array will be altered C if dry convective adjustment occurs and if IPBL is not equal to 0. C C QS: Array of saturation specific humidity of dimension ND, with first C index corresponding to lowest model level. Must be defined C at same grid levels as T. Note that this array will be altered C if dry convective adjustment occurs and if IPBL is not equal to 0. C C U: Array of zonal wind velocity (m/s) of dimension ND, witth first C index corresponding with the lowest model level. Defined at C same levels as T. Note that this array will be altered if C dry convective adjustment occurs and if IPBL is not equal to 0. C C V: Same as U but for meridional velocity. C C TRA: Array of passive tracer mixing ratio, of dimensions (ND,NTRA), C where NTRA is the number of different tracers. If no C convective tracer transport is needed, define a dummy C input array of dimension (ND,1). Tracers are defined at C same vertical levels as T. Note that this array will be altered C if dry convective adjustment occurs and if IPBL is not equal to 0. C C P: Array of pressure (mb) of dimension ND, with first C index corresponding to lowest model level. Must be defined C at same grid levels as T. C C PH: Array of pressure (mb) of dimension ND+1, with first index C corresponding to lowest level. These pressures are defined at C levels intermediate between those of P, T, Q and QS. The first C value of PH should be greater than (i.e. at a lower level than) C the first value of the array P. C C ND: The dimension of the arrays T,Q,QS,P,PH,FT and FQ C C NL: The maximum number of levels to which convection can C penetrate, plus 1. C NL MUST be less than or equal to ND-1. C C NTRA:The number of different tracers. If no tracer transport C is needed, set this equal to 1. (On most compilers, setting C NTRA to 0 will bypass tracer calculation, saving some CPU.) C C DELT: The model time step (sec) between calls to CONVECT C C -- sb: interface with the cloud parameterization: C C QCONDC: mixing ratio of condensed water within clouds (kg/kg) C For use in the Bony-Emanuel cloud parameterization C sb -- C C---------------------------------------------------------------------------- C *** On Output: *** C C IFLAG: An output integer whose value denotes the following: C C VALUE INTERPRETATION C ----- -------------- C 0 No moist convection; atmosphere is not C unstable, or surface temperature is less C than 250 K or surface specific humidity C is non-positive. C C 1 Moist convection occurs. C C 2 No moist convection: lifted condensation C level is above the 200 mb level. C C 3 No moist convection: cloud base is higher C then the level NL-1. C C 4 Moist convection occurs, but a CFL condition C on the subsidence warming is violated. This C does not cause the scheme to terminate. C C FT: Array of temperature tendency (K/s) of dimension ND, defined at same C grid levels as T, Q, QS and P. C C FQ: Array of specific humidity tendencies ((gm/gm)/s) of dimension ND, C defined at same grid levels as T, Q, QS and P. C C FU: Array of forcing of zonal velocity (m/s^2) of dimension ND, C defined at same grid levels as T. C C FV: Same as FU, but for forcing of meridional velocity. C C FTRA: Array of forcing of tracer content, in tracer mixing ratio per C second, defined at same levels as T. Dimensioned (ND,NTRA). C C PRECIP: Scalar convective precipitation rate (mm/day). C C WD: A convective downdraft velocity scale. For use in surface C flux parameterizations. See convect.ps file for details. C C TPRIME: A convective downdraft temperature perturbation scale (K). C For use in surface flux parameterizations. See convect.ps C file for details. C C QPRIME: A convective downdraft specific humidity C perturbation scale (gm/gm). C For use in surface flux parameterizations. See convect.ps C file for details. C C CBMF: The cloud base mass flux ((kg/m**2)/s). THIS SCALAR VALUE MUST C BE STORED BY THE CALLING PROGRAM AND RETURNED TO CONVECT AT C ITS NEXT CALL. That is, the value of CBMF must be "remembered" C by the calling program between calls to CONVECT. C C------------------------------------------------------------------------------ C C *** THE PARAMETER NA SHOULD IN GENERAL BE GREATER THAN *** C *** OR EQUAL TO ND + 1 *** C PARAMETER (NA=50) C INTEGER NENT(NA) REAL T(ND),Q(ND),QS(ND),U(ND),V(ND),TRA(ND,NTRA),P(ND),PH(ND) REAL FT(ND),FQ(ND),FU(ND),FV(ND),FTRA(ND,NTRA) REAL TRAENT(NA,NA,NTRA),TRATM(NA) REAL UP(NA),VP(NA),TRAP(NA,NTRA) REAL M(NA),MP(NA),MENT(NA,NA),QENT(NA,NA),ELIJ(NA,NA) REAL SIJ(NA,NA),TVP(NA),TV(NA),WATER(NA) REAL QP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA) REAL SIGP(NA),TP(NA),TOLD(NA),CPN(NA),LAMBDA REAL LV(NA),LVCP(NA),LV0,H(NA),HP(NA),GZ(NA),HM(NA) C -- ke 8/04 REAL UC(NA), VC(NA), UH(NA), VH(NA), FLUXU(NA), FLUXV(NA) C -- sb: REAL QCONDC(ND) REAL QCOND(NA),NQCOND(NA),WA(NA),MA(NA),SIGA(NA),AX(NA) C REAL MUP(NA),MDOWN(NA),ENT(NA),DET(NA),FQDET(NA),FQSUB(NA) COMMON / CVT / TVP,MP,QP,SIJ,MUP,MDOWN,ENT,DET,FQDET,FQSUB, 1 WATER, INB, CAPEM C C ----------------------------------------------------------------------- C C *** Specify Switches *** C C *** IPBL: Set to zero to bypass dry adiabatic adjustment *** C *** Any other value results in dry adiabatic adjustment *** C *** (Zero value recommended for use in models with *** C *** boundary layer schemes) *** C C *** MINORIG: Lowest level from which convection may originate *** C *** (Should be first model level at which T is defined *** C *** for models using bulk PBL schemes; otherwise, it should *** C *** be the first model level at which T is defined above *** C *** the surface layer) *** C IPBL=1 MINORIG=1 C C------------------------------------------------------------------------------ C C *** SPECIFY PARAMETERS *** C C *** ELCRIT IS THE AUTOCONVERSION THERSHOLD WATER CONTENT (gm/gm) *** C *** TLCRIT IS CRITICAL TEMPERATURE BELOW WHICH THE AUTO- *** C *** CONVERSION THRESHOLD IS ASSUMED TO BE ZERO *** C *** (THE AUTOCONVERSION THRESHOLD VARIES LINEARLY *** C *** BETWEEN 0 C AND TLCRIT) *** C *** ENTP IS THE COEFFICIENT OF MIXING IN THE ENTRAINMENT *** C *** FORMULATION *** C *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** C *** SIGS IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** C *** OF CLOUD *** C *** OMTRAIN IS THE ASSUMED FALL SPEED (P/s) OF RAIN *** C *** OMTSNOW IS THE ASSUMED FALL SPEED (P/s) OF SNOW *** C *** COEFFR IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** C *** OF RAIN *** C *** COEFFS IS A COEFFICIENT GOVERNING THE RATE OF EVAPORATION *** C *** OF SNOW *** C *** CU IS THE COEFFICIENT GOVERNING CONVECTIVE MOMENTUM *** C *** TRANSPORT *** C *** DTMAX IS THE MAXIMUM NEGATIVE TEMPERATURE PERTURBATION *** C *** A LIFTED PARCEL IS ALLOWED TO HAVE BELOW ITS LFC *** C *** ALPHA AND DAMP ARE PARAMETERS THAT CONTROL THE RATE OF *** C *** APPROACH TO QUASI-EQUILIBRIUM *** C *** (THEIR STANDARD VALUES ARE 0.20 AND 0.1, RESPECTIVELY) *** C *** (DAMP MUST BE LESS THAN 1) *** C ELCRIT=.0011 TLCRIT=-55.0 ENTP=1.5 SIGD=0.05 SIGS=0.12 SIGMIN=0.0 SIGMAX=0.9 OMTRAIN=50.0 OMTSNOW=5.5 COEFFR=0.9 COEFFS=0.6 LAMBDA=0.001/250.0 BETA=10.0 DTMAX=0.9 ALPHA=0.02 DAMP=0.1 DELTA=0.01 ! sb (for cloud parameterization) C c DTMAX=0.1*DTMAX1*(T(1)-273.15) c DTMAX=MIN(DTMAX,DTMAX1) c DTMAX=MAX(DTMAX,0.0) C C *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS, *** C *** GRAVITY, AND LIQUID WATER DENSITY. *** C *** THESE SHOULD BE CONSISTENT WITH *** C *** THOSE USED IN CALLING PROGRAM *** C *** NOTE: THESE ARE ALSO SPECIFIED IN SUBROUTINE TLIFT *** C CPD=1005.7 CPV=1870.0 CL=2500.0 RV=461.5 RD=287.04 LV0=2.501E6 G=9.8 ROWL=1000.0 C CPVMCL=CL-CPV EPS=RD/RV EPSI=1./EPS GINV=1.0/G DELTI=1.0/DELT C C *** INITIALIZE OUTPUT ARRAYS AND PARAMETERS *** C DO 5 I=1,ND FT(I)=0.0 FQ(I)=0.0 FU(I)=0.0 FV(I)=0.0 C -- sb: QCONDC(I)=0.0 QCOND(I)=0.0 NQCOND(I)=0.0 DO 4 J=1,NTRA FTRA(I,J)=0.0 4 CONTINUE MDOWN(I)=0.0 MUP(I)=0.0 ENT(I)=0.0 DET(I)=0.0 FQDET(I)=0.0 FQSUB(I)=0.0 WATER(I)=0.0 QP(I)=Q(I) MP(I)=0.0 TVP(I)=T(I)*(1.+Q(I)*EPSI-Q(I)) 5 CONTINUE DO 7 I=1,NL+1 RDCP=(RD*(1.-Q(I))+Q(I)*RV)/ 1 (CPD*(1.-Q(I))+Q(I)*CPV) TH(I)=T(I)*(1000.0/P(I))**RDCP 7 CONTINUE PRECIP=0.0 WD=0.0 TPRIME=0.0 QPRIME=0.0 IFLAG=0 INB=ND C C -- added by ke 8/04 C C *** INTERPOLATE TO FIND VELOCITIES AT HALF LEVELS *** C UH(ND)=U(ND) VH(ND)=V(ND) UC(ND)=UH(ND) VC(ND)=VH(ND) UP(ND)=UH(ND) VP(ND)=VH(ND) DO I=1,ND-1 UH(I)=0.5*(U(I)+U(I+1)) VH(I)=0.5*(V(I)+V(I+1)) UC(I)=UH(I) VC(I)=VH(I) UP(I)=UH(I) VP(I)=VH(I) END DO C C -- end of addition C IF(IPBL.NE.0)THEN C C *** PERFORM DRY ADIABATIC ADJUSTMENT *** C JC=0 DO 30 I=NL-1,1,-1 JN=0 SUM=TH(I)*(1.+Q(I)*EPSI-Q(I)) DO 10 J=I+1,NL SUM=SUM+TH(J)*(1.+Q(J)*EPSI-Q(J)) THBAR=SUM/FLOAT(J+1-I) IF((TH(J)*(1.+Q(J)*EPSI-Q(J))).LT.THBAR)JN=J IF(I.EQ.1.AND.P(J).GT.948.)JN=MAX(JN,J) 10 CONTINUE c IF(I.EQ.1)JN=MAX(JN,3) IF(JN.EQ.0)GOTO 30 12 CONTINUE AHM=0.0 RM=0.0 UM=0.0 VM=0.0 DO K=1,NTRA TRATM(K)=0.0 END DO DO 15 J=I,JN AHM=AHM+(CPD*(1.-Q(J))+Q(J)*CPV)*T(J)*(PH(J)-PH(J+1)) RM=RM+Q(J)*(PH(J)-PH(J+1)) UM=UM+U(J)*(PH(J)-PH(J+1)) VM=VM+V(J)*(PH(J)-PH(J+1)) DO K=1,NTRA TRATM(K)=TRATM(K)+TRA(J,K)*(PH(J)-PH(J+1)) END DO 15 CONTINUE DPHINV=1./(PH(I)-PH(JN+1)) RM=RM*DPHINV UM=UM*DPHINV VM=VM*DPHINV DO K=1,NTRA TRATM(K)=TRATM(K)*DPHINV END DO A2=0.0 DO 20 J=I,JN Q(J)=RM U(J)=UM V(J)=VM DO K=1,NTRA TRA(J,K)=TRATM(K) END DO RDCP=(RD*(1.-Q(J))+Q(J)*RV)/ 1 (CPD*(1.-Q(J))+Q(J)*CPV) X=(0.001*P(J))**RDCP TOLD(J)=T(J) T(J)=X A2=A2+(CPD*(1.-Q(J))+Q(J)*CPV)*X*(PH(J)-PH(J+1)) 20 CONTINUE DO 25 J=I,JN TH(J)=AHM/A2 T(J)=T(J)*TH(J) TC=TOLD(J)-273.15 ALV=LV0-CPVMCL*TC QS(J)=QS(J)+QS(J)*(1.+QS(J)*0.608)*ALV*(T(J)- 1 TOLD(J))/(RV*TOLD(J)*TOLD(J)) 25 CONTINUE IF((TH(JN+1)*(1.+Q(JN+1)*EPSI-Q(JN+1))).LT. 1 (TH(JN)*(1.+Q(JN)*EPSI-Q(JN))))THEN JN=JN+1 GOTO 12 END IF IF(I.EQ.1)JC=JN 30 CONTINUE C C *** Remove any supersaturation that results from adjustment *** C IF(JC.GT.1)THEN DO 38 J=1,JC IF(QS(J).LT.Q(J))THEN ALV=LV0-CPVMCL*(T(J)-273.15) TNEW=T(J)+ALV*(Q(J)-QS(J))/(CPD*(1.-Q(J))+ 1 CL*Q(J)+QS(J)*(CPV-CL+ALV*ALV/(RV*T(J)*T(J)))) ALVNEW=LV0-CPVMCL*(TNEW-273.15) QNEW=(ALV*Q(J)-(TNEW-T(J))*(CPD*(1.-Q(J))+CL*Q(J)))/ALVNEW PRECIP=PRECIP+24.*3600.*1.0E5*(PH(J)-PH(J+1))* 1 (Q(J)-QNEW)/(G*DELT*ROWL) T(J)=TNEW Q(J)=QNEW QS(J)=QNEW END IF 38 CONTINUE END IF C END IF C C *** CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY AND STATIC ENERGY C GZ(1)=0.0 CPN(1)=CPD*(1.-Q(1))+Q(1)*CPV H(1)=T(1)*CPN(1) LV(1)=LV0-CPVMCL*(T(1)-273.15) HM(1)=LV(1)*Q(1) TV(1)=T(1)*(1.+Q(1)*EPSI-Q(1)) AHMIN=1.0E12 IHMIN=NL DO 40 I=2,NL+1 TVX=T(I)*(1.+Q(I)*EPSI-Q(I)) TVY=T(I-1)*(1.+Q(I-1)*EPSI-Q(I-1)) GZ(I)=GZ(I-1)+0.5*RD*(TVX+TVY)*(P(I-1)-P(I))/PH(I) CPN(I)=CPD*(1.-Q(I))+CPV*Q(I) H(I)=T(I)*CPN(I)+GZ(I) LV(I)=LV0-CPVMCL*(T(I)-273.15) HM(I)=(CPD*(1.-Q(I))+CL*Q(I))*(T(I)-T(1))+ 1 LV(I)*Q(I)+GZ(I) TV(I)=T(I)*(1.+Q(I)*EPSI-Q(I)) C C *** Find level of minimum moist static energy *** C IF(I.GE.MINORIG.AND.HM(I).LT.AHMIN.AND.HM(I).LT.HM(I-1))THEN AHMIN=HM(I) IHMIN=I END IF 40 CONTINUE IHMIN=MIN(IHMIN, NL-1) C C *** Find that model level below the level of minimum moist *** C *** static energy that has the maximum value of moist static energy *** C AHMAX=0.0 DO 42 I=MINORIG,IHMIN IF(HM(I).GT.AHMAX)THEN NK=I AHMAX=HM(I) END IF 42 CONTINUE C C *** CHECK WHETHER PARCEL LEVEL TEMPERATURE AND SPECIFIC HUMIDITY *** C *** ARE REASONABLE *** C *** Skip convection if HM increases monotonically upward *** C IF(T(NK).LT.250.0.OR.Q(NK).LE.0.0.OR.IHMIN.EQ.(NL-1))THEN IFLAG=0 CBMF=0.0 RETURN END IF C C *** CALCULATE LIFTED CONDENSATION LEVEL OF AIR AT PARCEL ORIGIN LEVEL *** C *** (WITHIN 0.2% OF FORMULA OF BOLTON, MON. WEA. REV.,1980) *** C RH=Q(NK)/QS(NK) CHI=T(NK)/(1669.0-122.0*RH-T(NK)) PLCL=P(NK)*(RH**CHI) IF(PLCL.LT.200.0.OR.PLCL.GE.2000.0)THEN IFLAG=2 CBMF=0.0 RETURN END IF PBLH=(RD*T(NK)/G)*LOG(P(1)/PLCL) ROWS=100.*P(1)/(RD*T(1)) W3=PBLH*G*(FS/(CPD*T(1))+0.608*FL/2.5E6)/ROWS W3=MAX(W3,0.0) WSTAR=W3**(1./3.) DTMAX=1.737*WSTAR dtmax=max(dtmax,0.9) C C *** CALCULATE FIRST LEVEL ABOVE LCL (=ICB) *** C ICB=NL-1 DO 50 I=NK+1,NL IF(P(I).LT.PLCL)THEN ICB=MIN(ICB,I) END IF 50 CONTINUE IF(ICB.GE.(NL-1))THEN IFLAG=3 CBMF=0.0 RETURN END IF C C *** FIND TEMPERATURE UP THROUGH ICB AND TEST FOR INSTABILITY *** C C *** SUBROUTINE TLIFT CALCULATES PART OF THE LIFTED PARCEL VIRTUAL *** C *** TEMPERATURE, THE ACTUAL TEMPERATURE AND THE ADIABATIC *** C *** LIQUID WATER CONTENT *** C CALL TLIFT(P,T,Q,QS,GZ,ICB,NK,TVP,TP,CLW,ND,NL,1) DO 54 I=NK,ICB TVP(I)=TVP(I)-TP(I)*Q(NK) 54 CONTINUE C C *** If there was no convection at last time step and parcel *** C *** is stable at ICB then skip rest of calculation *** C IF(CBMF.EQ.0.0.AND.TVP(ICB).LE.(TV(ICB)-DTMAX))THEN IFLAG=0 RETURN END IF C C *** IF THIS POINT IS REACHED, MOIST CONVECTIVE ADJUSTMENT IS NECESSARY *** C IF(IFLAG.NE.4)IFLAG=1 C C *** FIND THE REST OF THE LIFTED PARCEL TEMPERATURES *** C CALL TLIFT(P,T,Q,QS,GZ,ICB,NK,TVP,TP,CLW,ND,NL,2) C C *** SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF *** C *** PRECIPITATION FALLING OUTSIDE OF CLOUD *** C *** THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) *** C DO 57 I=1,NK EP(I)=0.0 SIGP(I)=SIGS 57 CONTINUE DO 60 I=NK+1,NL TCA=TP(I)-273.15 IF(TCA.GE.0.0)THEN ELACRIT=ELCRIT ELSE ELACRIT=ELCRIT*(1.0-TCA/TLCRIT) END IF ELACRIT=MAX(ELACRIT,0.0) EPMAX=0.999 EP(I)=EPMAX*(1.0-ELACRIT/MAX(CLW(I),1.0E-8)) EP(I)=MAX(EP(I),0.0) EP(I)=MIN(EP(I),EPMAX) SIGP(I)=SIGS 60 CONTINUE C C *** CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL *** C *** VIRTUAL TEMPERATURE *** C DO 64 I=ICB+1,NL c TVP(I)=TVP(I)*(1.-Q(NK)+EP(I)*CLW(I)) TVP(I)=TVP(I)-Q(NK)*TP(I) 64 CONTINUE TVP(NL+1)=TVP(NL)-(GZ(NL+1)-GZ(NL))/CPD C C *** NOW INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS *** C DO 70 I=1,NL+1 HP(I)=H(I) NENT(I)=0 WATER(I)=0.0 EVAP(I)=0.0 WT(I)=OMTSNOW MP(I)=0.0 M(I)=0.0 LVCP(I)=LV(I)/CPN(I) C -- ke 8/04 FLUXU(I)=0.0 FLUXV(I)=0.0 C -- DO 70 J=1,NL+1 QENT(I,J)=Q(J) ELIJ(I,J)=0.0 MENT(I,J)=0.0 SIJ(I,J)=0.0 DO 70 K=1,NTRA TRAENT(I,J,K)=TRA(J,K) 70 CONTINUE QP(1)=Q(1) DO 71 I=1,NTRA TRAP(1,I)=TRA(1,I) 71 CONTINUE DO 72 I=2,NL+1 QP(I)=Q(I-1) DO 72 J=1,NTRA TRAP(I,J)=TRA(I-1,J) 72 CONTINUE C C *** FIND THE FIRST MODEL LEVEL (INB1) ABOVE THE PARCEL'S *** C *** HIGHEST LEVEL OF NEUTRAL BUOYANCY *** C *** AND THE HIGHEST LEVEL OF POSITIVE CAPE (INB) *** C CAPE=0.0 CAPEM=0.0 INB=ICB+1 INB1=INB BYP=0.0 DO 82 I=ICB+1,NL-1 BY=(TVP(I)-TV(I))*(PH(I)-PH(I+1))/P(I) CAPE=CAPE+BY IF(BY.GE.0.0)INB1=I+1 IF(CAPE.GT.0.0)THEN INB=I+1 CAPEM=CAPE BYP=(TVP(I+1)-TV(I+1))*(PH(I+1)-PH(I+2))/P(I+1) END IF 82 CONTINUE c INB=MAX(INB,INB1) INB=(INB+INB1)/2 CAPE=CAPEM+BYP DEFRAC=CAPEM-CAPE DEFRAC=MAX(DEFRAC,0.001) FRAC=-CAPE/DEFRAC FRAC=MIN(FRAC,1.0) FRAC=MAX(FRAC,0.0) JMAX=INB C C *** CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL *** C DO 95 I=ICB,INB HP(I)=H(NK)+(LV(I)+(CPD-CPV)*T(I))*EP(I)*CLW(I) 95 CONTINUE C C *** CALCULATE CLOUD BASE MASS FLUX AND RATES OF MIXING, M(I), *** c *** AT EACH MODEL LEVEL *** C DBOSUM=0.0 C C *** INTERPOLATE DIFFERENCE BETWEEN LIFTED PARCEL AND *** C *** ENVIRONMENTAL TEMPERATURES TO LIFTED CONDENSATION LEVEL *** C TVPPLCL=TVP(ICB-1)-RD*TVP(ICB-1)*(P(ICB-1)-PLCL)/ 1 (CPN(ICB-1)*P(ICB-1)) TVAPLCL=TV(ICB)+(TVP(ICB)-TVP(ICB+1))*(PLCL-P(ICB))/ 1 (P(ICB)-P(ICB+1)) DTPBL=0.0 DO 96 I=NK,ICB-1 DTPBL=DTPBL+(TVP(I)-TV(I))*(PH(I)-PH(I+1)) 96 CONTINUE DTPBL=DTPBL/(PH(NK)-PH(ICB)) DTMIN=TVPPLCL-TVAPLCL+DTMAX+DTPBL DTMA=DTMIN C C *** ADJUST CLOUD BASE MASS FLUX *** C CBMFOLD=CBMF DELT0=300.0 DAMPS=DAMP*DELT/DELT0 CBMF=(1.-DAMPS)*CBMF+0.1*ALPHA*DTMA c CBMF=0.2*ALPHA*DTMA/DAMPS CBMF=MAX(CBMF,0.0) C C *** If cloud base mass flux is zero, skip rest of calculation *** C IF(CBMF.EQ.0.0.AND.CBMFOLD.EQ.0.0)THEN RETURN END IF C C *** CALCULATE RATES OF MIXING, M(I) *** C M(ICB)=0.0 DO 103 I=ICB+1,INB K=MIN(I,INB1) c DBO=ABS(TV(K+1)-TVP(K+1)-TV(K-1)+TVP(K-1))+ c 1 ENTP*0.04*(PH(K)-PH(K+1)) c c -- modif proposed at NRL (Monterey): DBO=ABS(TV(K)-TVP(K))+ 1 ENTP*0.02*(PH(K)-PH(K+1)) c DBOSUM=DBOSUM+DBO M(I)=CBMF*DBO 103 CONTINUE DO 110 I=ICB+1,INB M(I)=M(I)/DBOSUM 110 CONTINUE C C *** CALCULATE ENTRAINED AIR MASS FLUX (MENT), TOTAL WATER MIXING *** C *** RATIO (QENT), TOTAL CONDENSED WATER (ELIJ), AND MIXING *** C *** FRACTION (SIJ) *** C DO 170 I=ICB+1,INB QTI=Q(NK)-EP(I)*CLW(I) DO 160 J=ICB,JMAX BF2=1.+LV(J)*LV(J)*QS(J)/(RV*T(J)*T(J)*CPD) ANUM=H(J)-HP(I)+(CPV-CPD)*T(J)*(QTI-Q(J)) c ANUM=H(J)-HP(I)+(CPV-CPD)*T(J)*(QTI-Q(J))-CPN(J)*Q(NK)*TP(J) DENOM=H(I)-HP(I)+(CPD-CPV)*(Q(I)-QTI)*T(J) DEI=DENOM IF(ABS(DEI).LT.0.01)DEI=0.01 SIJ(I,J)=ANUM/DEI SIJ(I,I)=1.0 ALTEM=SIJ(I,J)*Q(I)+(1.-SIJ(I,J))*QTI-QS(J) ALTEM=ALTEM/BF2 CWAT=CLW(J)*(1.-EP(J)) STEMP=SIJ(I,J) IF((STEMP.LT.0.0.OR.STEMP.GT.1.0.OR. 1 ALTEM.GT.CWAT).AND.J.GT.I)THEN ANUM=ANUM-LV(J)*(QTI-QS(J)-CWAT*BF2) DENOM=DENOM+LV(J)*(Q(I)-QTI) IF(ABS(DENOM).LT.0.01)DENOM=0.01 SIJ(I,J)=ANUM/DENOM ALTEM=SIJ(I,J)*Q(I)+(1.-SIJ(I,J))*QTI-QS(J) ALTEM=ALTEM-(BF2-1.)*CWAT END IF IF(SIJ(I,J).GT.SIGMIN.AND.SIJ(I,J).LT.SIGMAX)THEN QENT(I,J)=SIJ(I,J)*Q(I)+(1.-SIJ(I,J))*QTI DO K=1,NTRA TRAENT(I,J,K)=SIJ(I,J)*TRA(I,K)+(1.-SIJ(I,J))* 1 TRA(NK,K) END DO ELIJ(I,J)=ALTEM ELIJ(I,J)=MAX(0.0,ELIJ(I,J)) MENT(I,J)=M(I)/(1.-SIJ(I,J)) NENT(I)=NENT(I)+1 END IF SIJ(I,J)=MAX(0.0,SIJ(I,J)) SIJ(I,J)=MIN(1.0,SIJ(I,J)) 160 CONTINUE C C *** IF NO AIR CAN ENTRAIN AT LEVEL I ASSUME THAT UPDRAFT DETRAINS *** C *** AT THAT LEVEL AND CALCULATE DETRAINED AIR FLUX AND PROPERTIES *** C IF(NENT(I).EQ.0)THEN MENT(I,I)=M(I) QENT(I,I)=Q(NK)-EP(I)*CLW(I) DO J=1,NTRA TRAENT(I,I,J)=TRA(NK,J) END DO ELIJ(I,I)=CLW(I) SIJ(I,I)=1.0 END IF 170 CONTINUE SIJ(INB,INB)=1.0 C C *** NORMALIZE ENTRAINED AIR MASS FLUXES TO REPRESENT EQUAL *** C *** PROBABILITIES OF MIXING *** C DO 200 I=ICB+1,INB IF(NENT(I).NE.0)THEN QP1=Q(NK)-EP(I)*CLW(I) ANUM=H(I)-HP(I)-LV(I)*(QP1-QS(I)) DENOM=H(I)-HP(I)+LV(I)*(Q(I)-QP1) IF(ABS(DENOM).LT.0.01)DENOM=0.01 SCRIT=ANUM/DENOM ALT=QP1-QS(I)+SCRIT*(Q(I)-QP1) IF(ALT.LT.0.0)SCRIT=1.0 SCRIT=MAX(SCRIT,0.0) ASIJ=0.0 SMIN=1.0 DO 175 J=ICB,JMAX IF(SIJ(I,J).GT.SIGMIN.AND.SIJ(I,J).LT.SIGMAX)THEN IF(J.GT.I)THEN SMID=MIN(SIJ(I,J),SCRIT) SJMAX=SMID SJMIN=SMID IF(SMID.LT.SMIN.AND.SIJ(I,J+1).LT.SMID)THEN SMIN=SMID SJMAX=MIN(SIJ(I,J+1),SIJ(I,J),SCRIT) SJMIN=MAX(SIJ(I,J-1),SIJ(I,J)) SJMIN=MIN(SJMIN,SCRIT) END IF ELSE SJMAX=MAX(SIJ(I,J+1),SCRIT) SMID=MAX(SIJ(I,J),SCRIT) SJMIN=0.0 IF(J.GT.1)SJMIN=SIJ(I,J-1) SJMIN=MAX(SJMIN,SCRIT) END IF DELP=ABS(SJMAX-SMID) DELM=ABS(SJMIN-SMID) ASIJ=ASIJ+(DELP+DELM)*(PH(J)-PH(J+1)) MENT(I,J)=MENT(I,J)*(DELP+DELM)*(PH(J)-PH(J+1)) END IF 175 CONTINUE ASIJ=MAX(1.0E-21,ASIJ) ASIJ=1.0/ASIJ DO 180 J=ICB,INB MENT(I,J)=MENT(I,J)*ASIJ 180 CONTINUE BSUM=0.0 DO 190 J=ICB,INB BSUM=BSUM+MENT(I,J) 190 CONTINUE IF(BSUM.LT.1.0E-18)THEN NENT(I)=0 MENT(I,I)=M(I) QENT(I,I)=Q(NK)-EP(I)*CLW(I) DO J=1,NTRA TRAENT(I,I,J)=TRA(NK,J) END DO ELIJ(I,I)=CLW(I) SIJ(I,I)=1.0 END IF END IF 200 CONTINUE C C *** CHECK WHETHER EP(INB)=0, IF SO, SKIP PRECIPITATING *** C *** DOWNDRAFT CALCULATION *** C IF(EP(INB).LT.0.0001)GOTO 405 C C *** INTEGRATE LIQUID WATER EQUATION TO FIND CONDENSED WATER *** C *** AND CONDENSED WATER FLUX *** C JTT=2 C C *** BEGIN DOWNDRAFT LOOP *** C DO 400 I=INB,1,-1 C C *** CALCULATE DETRAINED PRECIPITATION *** C WDTRAIN=G*EP(I)*M(I)*CLW(I) IF(I.GT.1)THEN DO 320 J=1,I-1 AWAT=ELIJ(J,I)-(1.-EP(I))*CLW(I) AWAT=MAX(0.0,AWAT) 320 WDTRAIN=WDTRAIN+G*AWAT*MENT(J,I) END IF C C *** FIND RAIN WATER AND EVAPORATION USING PROVISIONAL *** C *** ESTIMATES OF QP(I)AND QP(I-1) *** C c c *** Value of terminal velocity and coeffecient of evaporation for snow *** c COEFF=COEFFS WT(I)=OMTSNOW c c *** Value of terminal velocity and coeffecient of evaporation for rain *** c IF(T(I).GT.273.0)THEN COEFF=COEFFR WT(I)=OMTRAIN END IF QSM=0.5*(Q(I)+QP(I+1)) AFAC=COEFF*PH(I)*(QS(I)-QSM)/(1.0E4+2.0E3*PH(I)*QS(I)) AFAC=MAX(AFAC,0.0) SIGT=SIGP(I) SIGT=MAX(0.0,SIGT) SIGT=MIN(1.0,SIGT) B6=100.*(PH(I)-PH(I+1))*SIGT*AFAC/WT(I) C6=(WATER(I+1)*WT(I+1)+WDTRAIN/SIGD)/WT(I) REVAP=0.5*(-B6+SQRT(B6*B6+4.*C6)) EVAP(I)=SIGT*AFAC*REVAP WATER(I)=REVAP*REVAP C C *** CALCULATE PRECIPITATING DOWNDRAFT MASS FLUX UNDER *** C *** HYDROSTATIC APPROXIMATION *** C IF(I.EQ.1)GOTO 360 DHDP=(H(I)-H(I-1))/(P(I-1)-P(I)) DHDP=MAX(DHDP,10.0) MP(I)=100.*GINV*LV(I)*SIGD*EVAP(I)/DHDP MP(I)=MAX(MP(I),0.0) C C *** ADD SMALL AMOUNT OF INERTIA TO DOWNDRAFT *** C FAC=20.0/(PH(I-1)-PH(I)) MP(I)=(FAC*MP(I+1)+MP(I))/(1.+FAC) C C *** FORCE MP TO DECREASE LINEARLY TO ZERO *** C *** BETWEEN ABOUT 900 MB AND THE SURFACE *** C IF(P(I).GT.(0.899*P(1)))THEN JTT=MAX(JTT,I) MP(I)=MP(JTT)*(P(1)-P(I))/(P(1)-P(JTT)) END IF 360 CONTINUE C C *** FIND MIXING RATIO OF PRECIPITATING DOWNDRAFT *** C IF(I.EQ.INB)GOTO 400 IF(I.EQ.1)THEN QSTM=QS(1) ELSE QSTM=QS(I-1) END IF IF(MP(I).GT.MP(I+1))THEN RAT=MP(I+1)/MP(I) QP(I)=QP(I+1)*RAT+Q(I)*(1.0-RAT)+100.*GINV* 1 SIGD*(PH(I)-PH(I+1))*(EVAP(I)/MP(I)) DO J=1,NTRA TRAP(I,J)=TRAP(I+1,J)*RAT+TRAP(I,J)*(1.-RAT) END DO ELSE IF(MP(I+1).GT.0.0)THEN QP(I)=(GZ(I+1)-GZ(I)+QP(I+1)*(LV(I+1)+T(I+1)*( 1 CL-CPD))+CPD*(T(I+1)-T(I)))/(LV(I)+T(I)*(CL-CPD)) DO J=1,NTRA TRAP(I,J)=TRAP(I+1,J) END DO END IF END IF QP(I)=MIN(QP(I),QSTM) QP(I)=MAX(QP(I),0.0) 400 CONTINUE C C *** CALCULATE SURFACE PRECIPITATION IN MM/DAY *** C PRECIP=PRECIP+WT(1)*SIGD*WATER(1)*3600.*24000./(ROWL*G) C 405 CONTINUE C C *** CALCULATE DOWNDRAFT VELOCITY SCALE AND SURFACE TEMPERATURE AND *** c *** WATER VAPOR FLUCTUATIONS *** C WD=BETA*ABS(MP(ICB))*0.01*RD*T(ICB)/(SIGD*P(ICB)) QPRIME=0.5*(QP(1)-Q(1)) TPRIME=LV0*QPRIME/CPD C C -- ke 8/04 C C *** Calculate in-cloud horizontal velocities and fluxes *** C ICBM=MAX(ICB,3) C IF(INB.GE.ICBM)THEN C DO I=ICBM,INB AMP1=0.0 IF(I.GE.NK)THEN DO K=I,INB AMP1=AMP1+M(K) END DO END IF DO K=1,I DO J=I,INB AMP1=AMP1+MENT(K,J) END DO END DO DO K=1,I-1 DO J=I,INB AMP1=AMP1-MENT(J,K) END DO END DO AMPN=MAX(AMP1,0.0) AMP1=MAX(AMP1,1.0E-4) C CHIZ=50.*(PH(I-1)-PH(I+1))*LAMBDA/(AMP1*G) CHIU=CHIZ*ABS(UC(I-1)-UH(I-1)) UC(I)=(UC(I-2)*(1.-CHIU)+CHIU*(UH(I)+UH(I-2)))/(1.+CHIU) UC(I-1)=UC(I-1)+0.4*(UC(I)+UC(I-2)-2.*UC(I-1)) FLUXU(I-1)=AMPN*(UC(I-1)-UH(I-1)) CHIV=CHIZ*ABS(VC(I-1)-VH(I-1)) VC(I)=(VC(I-2)*(1.-CHIV)+CHIV*(VH(I)+VH(I-2)))/(1.+CHIV) VC(I-1)=VC(I-1)+0.4*(VC(I)+VC(I-2)-2.*VC(I-1)) FLUXV(I-1)=AMPN*(VC(I-1)-VH(I-1)) END DO IF(ICB.GT.(NK+1))THEN DO I=NK+1,ICB-1 DFAC=FLOAT(I-NK)/FLOAT(ICB-NK) FLUXU(I)=FLUXU(ICB)*DFAC FLUXV(I)=FLUXV(ICB)*DFAC END DO END IF C END IF C C *** Calculate unsaturated downdraft momentum flux *** C IF(INB.GT.2)THEN C DO I=INB-2,1,-1 IF(MP(I+2).GT.0.0)THEN AMP1=MAX(MP(I+2),1.0E-4) CHIZ=50.*(PH(I+1)-PH(I+3))*LAMBDA/(G*AMP1) CHIU=CHIZ*ABS(UP(I+1)-UH(I+1)) UP(I)=(UP(I+2)*(1.-CHIU)+CHIU*(UH(I)+UH(I+2)))/(1.+CHIU) UP(I+1)=UP(I+1)+0.4*(UP(I)+UP(I+2)-2.*UP(I+1)) FLUXU(I+1)=FLUXU(I+1)-MP(I+2)*(UP(I+1)-UH(I+1)) CHIV=CHIZ*ABS(VP(I+1)-VH(I+1)) VP(I)=(VP(I+2)*(1.-CHIV)+CHIV*(VH(I)+VH(I+2)))/(1.+CHIV) VP(I+1)=VP(I+1)+0.4*(VP(I)+VP(I+2)-2.*VP(I+1)) FLUXV(I+1)=FLUXV(I+1)-MP(I+2)*(VP(I+1)-VH(I+1)) END IF END DO FLUXU(1)=-MP(2)*(UP(1)-UH(1)) FLUXV(1)=-MP(2)*(VP(1)-VH(1)) C END IF C C -- end of addition by ke 8/04 C C *** CALCULATE TENDENCIES OF LOWEST LEVEL POTENTIAL TEMPERATURE *** C *** AND MIXING RATIO *** C DPINV=0.01/(PH(1)-PH(2)) AM=0.0 IF(NK.EQ.1)THEN DO 410 K=2,INB 410 AM=AM+M(K) END IF IF((2.*G*DPINV*AM).GE.DELTI)IFLAG=4 FT(1)=FT(1)+G*DPINV*AM*(T(2)-T(1)+(GZ(2)-GZ(1))/CPN(1)) FT(1)=FT(1)-LVCP(1)*SIGD*EVAP(1) FT(1)=FT(1)+SIGD*WT(2)*(CL-CPD)*WATER(2)*(T(2)- 1 T(1))*DPINV/CPN(1) c FT(1)=FT(1)+0.01*SIGD*WT(1)*WATER(1)*RD*TV(1)/(PH(1)*CPN(1)) FQ(1)=FQ(1)+G*MP(2)*(QP(2)-Q(1))* 1 DPINV+SIGD*EVAP(1) FQDET(1)=FQ(1) FQ(1)=FQ(1)+G*AM*(Q(2)-Q(1))*DPINV FQSUB(1)=G*AM*(Q(2)-Q(1))*DPINV FU(1)=-G*DPINV*FLUXU(1) FV(1)=-G*DPINV*FLUXV(1) DO J=1,NTRA FTRA(1,J)=FTRA(1,J)+G*DPINV*(MP(2)*(TRAP(2,J)-TRA(1,J))+ 1 AM*(TRA(2,J)-TRA(1,J))) END DO AMDE=0.0 DO 415 J=2,INB FQ(1)=FQ(1)+G*DPINV*MENT(J,1)*(QENT(J,1)-Q(1)) FQSUB(1)=FQSUB(1)+G*DPINV*MENT(J,1)*(QENT(J,1)-Q(1)) DO K=1,NTRA FTRA(1,K)=FTRA(1,K)+G*DPINV*MENT(J,1)*(TRAENT(J,1,K)- 1 TRA(1,K)) END DO 415 CONTINUE C C *** CALCULATE TENDENCIES OF POTENTIAL TEMPERATURE AND MIXING RATIO *** C *** AT LEVELS ABOVE THE LOWEST LEVEL *** C C *** FIRST FIND THE NET SATURATED UPDRAFT AND DOWNDRAFT MASS FLUXES *** C *** THROUGH EACH LEVEL *** C DO 500 I=2,INB DPINV=0.01/(PH(I)-PH(I+1)) CPINV=1.0/CPN(I) AMP1=0.0 AD=0.0 IF(I.GE.NK)THEN c IF(I.GE.ICB-1)THEN DO 440 K=I+1,INB+1 440 AMP1=AMP1+M(K) END IF DO 450 K=1,I DO 450 J=I+1,INB+1 AMP1=AMP1+MENT(K,J) IF(I.GE.ICB)MUP(I)=AMP1 c MUP(I)=AMP1 450 CONTINUE IF((2.*G*DPINV*AMP1).GE.DELTI)THEN IFLAG=4 c CBMF=CBMFOLD-0.1*(2.*G*DPINV*AMP1*DELT-1.) c CBMF=MAX(CBMF,0.0) END IF DO 470 K=1,I-1 DO 470 J=I,INB AD=AD+MENT(J,K) MDOWN(I)=AD 470 CONTINUE FT(I)=FT(I)+G*DPINV*(AMP1*(T(I+1)-T(I)+(GZ(I+1)-GZ(I))* 1 CPINV)-AD*(T(I)-T(I-1)+(GZ(I)-GZ(I-1))*CPINV)) 2 -SIGD*LVCP(I)*EVAP(I) FT(I)=FT(I)+G*DPINV*MENT(I,I)*(HP(I)-H(I)+ 1 T(I)*(CPV-CPD)*(Q(I)-QENT(I,I)))*CPINV FT(I)=FT(I)+SIGD*WT(I+1)*(CL-CPD)*WATER(I+1)* 1 (T(I+1)-T(I))*DPINV*CPINV FQ(I)=FQ(I)+G*DPINV*(AMP1*(Q(I+1)-Q(I))- 1 AD*(Q(I)-Q(I-1))) C C -- added by ke 8/04 C C *** Momentum fluxes *** C FU(I)=-G*DPINV*(FLUXU(I)-FLUXU(I-1)) FV(I)=-G*DPINV*(FLUXV(I)-FLUXV(I-1)) C C -- end of addition C DO K=1,NTRA FTRA(I,K)=FTRA(I,K)+G*DPINV*(AMP1*(TRA(I+1,K)- 1 TRA(I,K))-AD*(TRA(I,K)-TRA(I-1,K))) END DO FQSUB(I)=FQ(I) DO 480 K=1,I-1 AWAT=ELIJ(K,I)-(1.-EP(I))*CLW(I) AWAT=MAX(AWAT,0.0) FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-AWAT-Q(I)) FQDET(I)=G*DPINV*MENT(K,I)*(QENT(K,I)-AWAT-Q(I)) c -- sb: C (saturated updrafts resulting from mixing) QCOND(I)=QCOND(I)+(ELIJ(K,I)-AWAT) NQCOND(I)=NQCOND(I)+1. c sb -- DO J=1,NTRA FTRA(I,J)=FTRA(I,J)+G*DPINV*MENT(K,I)*(TRAENT(K,I,J)- 1 TRA(I,J)) END DO 480 CONTINUE DO 490 K=I,INB FQ(I)=FQ(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-Q(I)) FQDET(I)=FQDET(I)+G*DPINV*MENT(K,I)*(QENT(K,I)-Q(I)) DO J=1,NTRA FTRA(I,J)=FTRA(I,J)+G*DPINV*MENT(K,I)*(TRAENT(K,I,J)- 1 TRA(I,J)) END DO 490 CONTINUE FQ(I)=FQ(I)+SIGD*EVAP(I)+G*(MP(I+1)* 1 (QP(I+1)-Q(I))-MP(I)*(QP(I)-Q(I-1)))*DPINV FQDET(I)=FQDET(I)+SIGD*EVAP(I)+G*(MP(I+1)* 1 (QP(I+1)-Q(I))-MP(I)*(QP(I)-Q(I-1)))*DPINV DO J=1,NTRA FTRA(I,J)=FTRA(I,J)+G*DPINV*(MP(I+1)*(TRAP(I+1,J)-TRA(I,J))- 1 MP(I)*(TRAP(I,J)-TRAP(I-1,J))) END DO c -- sb: C (saturated downdrafts resulting from mixing) DO K=I+1,INB QCOND(I)=QCOND(I)+ELIJ(K,I) NQCOND(I)=NQCOND(I)+1. ENDDO C (particular case: no detraining level is found) IF (NENT(I).EQ.0) THEN QCOND(I)=QCOND(I)+(1-EP(I))*CLW(I) NQCOND(I)=NQCOND(I)+1. ENDIF IF (NQCOND(I).NE.0.) THEN QCOND(I)=QCOND(I)/NQCOND(I) ENDIF c sb -- 500 CONTINUE C C *** Adjust tendencies at top of convection layer to reflect *** C *** actual position of the level zero CAPE *** C FQOLD=FQ(INB) FQ(INB)=FQ(INB)*(1.-FRAC) FQ(INB-1)=FQ(INB-1)+FRAC*FQOLD*((PH(INB)-PH(INB+1))/ 1 (PH(INB-1)-PH(INB)))*LV(INB)/LV(INB-1) FTOLD=FT(INB) FT(INB)=FT(INB)*(1.-FRAC) FT(INB-1)=FT(INB-1)+FRAC*FTOLD*((PH(INB)-PH(INB+1))/ 1 (PH(INB-1)-PH(INB)))*CPN(INB)/CPN(INB-1) FUOLD=FU(INB) FU(INB)=FU(INB)*(1.-FRAC) FU(INB-1)=FU(INB-1)+FRAC*FUOLD*((PH(INB)-PH(INB+1))/ 1 (PH(INB-1)-PH(INB))) FVOLD=FV(INB) FV(INB)=FV(INB)*(1.-FRAC) FV(INB-1)=FV(INB-1)+FRAC*FVOLD*((PH(INB)-PH(INB+1))/ 1 (PH(INB-1)-PH(INB))) DO K=1,NTRA FTRAOLD=FTRA(INB,K) FTRA(INB,K)=FTRA(INB,K)*(1.-FRAC) FTRA(INB-1,K)=FTRA(INB-1,K)+FRAC*FTRAOLD*(PH(INB)-PH(INB+1))/ 1 (PH(INB-1)-PH(INB)) END DO C C *** Very slightly adjust tendencies to force exact *** C *** enthalpy, momentum and tracer conservation *** C ENTS=0.0 UAV=0.0 VAV=0.0 DO 680 I=1,INB ENTS=ENTS+(CPN(I)*FT(I)+LV(I)*FQ(I))*(PH(I)-PH(I+1)) UAV=UAV+FU(I)*(PH(I)-PH(I+1)) VAV=VAV+FV(I)*(PH(I)-PH(I+1)) 680 CONTINUE ENTS=ENTS/(PH(1)-PH(INB+1)) UAV=UAV/(PH(1)-PH(INB+1)) VAV=VAV/(PH(1)-PH(INB+1)) DO 640 I=1,INB FT(I)=FT(I)-ENTS/CPN(I) FU(I)=FU(I)-UAV FV(I)=FV(I)-VAV 640 CONTINUE DO 700 K=1,NTRA TRAAV=0.0 DO 690 I=1,INB TRAAV=TRAAV+FTRA(I,K)*(PH(I)-PH(I+1)) 690 CONTINUE TRAAV=TRAAV/(PH(1)-PH(INB+1)) DO 695 I=1,INB FTRA(I,K)=FTRA(I,K)-TRAAV 695 CONTINUE 700 CONTINUE DO 750 I=1,INB DO 710 K=1,INB DET(I)=DET(I)+MENT(K,I) 710 CONTINUE DO 720 K=1,INB ENT(I)=ENT(I)+MENT(I,K) 720 CONTINUE ENT(I)=ENT(I)-MENT(I,I) 750 CONTINUE C In-cloud mixing ratio of condensed water : DO I=1,ND MA(I)=0.0 WA(I)=0.0 SIGA(I)=0.0 ENDDO DO I=NK,INB DO K=I+1,INB+1 MA(I)=MA(I)+M(K) ENDDO ENDDO DO I=ICB,INB-1 AX(I)=0. DO J=ICB,I AX(I)=AX(I)+RD*(TVP(J)-TV(J))*(PH(J)-PH(J+1))/P(J) ENDDO IF (AX(I).GT.0.) THEN WA(I)=SQRT(2.*AX(I)) ENDIF ENDDO DO I=1,NL IF (WA(I).GT.0. ) : SIGA(I)=MA(I)/WA(I)*RD*TVP(I)/P(I)/100./DELTA SIGA(I) = MIN(SIGA(I),1.0) QCONDC(I)=SIGA(I)*CLW(I)*(1.-EP(I)) : + (1.-SIGA(I))*QCOND(I) ENDDO C C *** RETURN *** C RETURN C END C C --------------------------------------------------------------------------- C SUBROUTINE TLIFT(P,T,Q,QS,GZ,ICB,NK,TVP,TPK,CLW,ND,NL,KK) REAL GZ(ND),TPK(ND),CLW(ND),P(ND) REAL T(ND),Q(ND),QS(ND),TVP(ND),LV0 C C *** ASSIGN VALUES OF THERMODYNAMIC CONSTANTS *** C CPD=1005.7 CPV=1870.0 CL=2500.0 RV=461.5 RD=287.04 LV0=2.501E6 C CPVMCL=CL-CPV EPS=RD/RV EPSI=1./EPS C C *** CALCULATE CERTAIN PARCEL QUANTITIES, INCLUDING STATIC ENERGY *** C ALV=LV0-CPVMCL*(T(NK)-273.15) AH0=(CPD*(1.-Q(NK))+CL*Q(NK))*T(NK)+Q(NK)*ALV+GZ(NK) CPP=CPD*(1.-Q(NK))+Q(NK)*CPV CPINV=1./CPP C IF(KK.EQ.1)THEN C C *** CALCULATE LIFTED PARCEL QUANTITIES BELOW CLOUD BASE *** C DO 50 I=1,ICB-1 CLW(I)=0.0 50 CONTINUE DO 100 I=NK,ICB-1 TPK(I)=T(NK)-(GZ(I)-GZ(NK))*CPINV TVP(I)=TPK(I)*(1.+Q(NK)*EPSI) 100 CONTINUE END IF C C *** FIND LIFTED PARCEL QUANTITIES ABOVE CLOUD BASE *** C NST=ICB NSB=ICB IF(KK.EQ.2)THEN NST=NL NSB=ICB+1 END IF DO 300 I=NSB,NST TG=T(I) QG=QS(I) ALV=LV0-CPVMCL*(T(I)-273.15) DO 200 J=1,2 S=CPD+ALV*ALV*QG/(RV*T(I)*T(I)) S=1./S AHG=CPD*TG+(CL-CPD)*Q(NK)*T(I)+ALV*QG+GZ(I) TG=TG+S*(AH0-AHG) TG=MAX(TG,35.0) TC=TG-273.15 DENOM=243.5+TC IF(TC.GE.0.0)THEN ES=6.112*EXP(17.67*TC/DENOM) ELSE ES=EXP(23.33086-6111.72784/TG+0.15215*LOG(TG)) END IF QG=EPS*ES/(P(I)-ES*(1.-EPS)) 200 CONTINUE ALV=LV0-CPVMCL*(T(I)-273.15) TPK(I)=(AH0-(CL-CPD)*Q(NK)*T(I)-GZ(I)-ALV*QG)/CPD CLW(I)=Q(NK)-QG CLW(I)=MAX(0.0,CLW(I)) RG=QG/(1.-Q(NK)) TVP(I)=TPK(I)*(1.+RG*EPSI) 300 CONTINUE RETURN END C SUBROUTINE CLOUDS_SUB_LS_40(ND,R,RS,T,P,PH,DT,QSUBGRID : ,CLDF,CLDQ,PRADJ,FTADJ,FRADJ,SIGSUB) implicit none C C========================================================================== C C CLOUDS_SUB_LS version 4.0 C C Purpose: C -------- C C 1) Call of the cloud parameterization C C 2) Large-scale super-saturation adjustment: C - condense water that exceeds saturation C - precipitate a fraction of that condensed water C - adjust the final in-cloud water content C C Inputs: C ------- C ND----------: Number of vertical levels C R--------ND-: Grid-box average of the total water mixing ratio [kg/kg] C RS-------ND-: Mean saturation humidity mixing ratio within the gridbox [kg/kg] C T--------ND-: Grid-box average temperature [K] C P-------ND+1: Pressure at mid-levels [mb] C PH------ND+1: Pressure at interface levels [mb] C DT----------: Timestep [seconds] C QSUBGRID-ND-: in-cloud mixing ratio of cloud condensate [kg/kg] from CONVECT C C Outputs: C -------- C CLDF-----ND-: cloud fraction [0-1] C CLDQ-----ND-: in-cloud mixing ratio of condensed water [kg/kg] C PRADJ----ND-: precipitation associated with the LS super-saturation [mm/day] C FTADJ----ND-: temperature tendency associated with the LS adjustment [K/s] C FRADJ----ND-: total water tendency associated with the LS adjustment [kg/kg/s] C SIGSUB---ND-: SQRT(variance) of total water [kg/kg] (diagnostic only) C C Written by: C ----------- C Sandrine Bony (MIT & LMD/CNRS) - August 1999 C C Modified by: C ------------ C C============================================================================ integer NA,ND,I parameter (NA=50) real DT, PRADJ real TCA,ELACRIT,ALV,CPN,TNEW,RNEW,EP real R(ND),RS(ND),T(ND),P(ND),PH(ND+1) : ,QSUBGRID(ND),CLDF(ND),CLDQ(ND) : ,FTADJ(ND),FRADJ(ND),RNEWLS(NA),TNEWLS(NA),QLSP(NA) : ,EPLS(NA) c a retirer dans version pub: real RSSURF,CLDL(NA),CLDA(NA),CLDK(NA),CLDS(NA) real BIDON(ND) real SIGSUB(ND) COMMON / cldpara / CLDL, CLDA, CLDK, CLDS C-------------------------------------------------------------------- C Thermodynamical constants: C REAL CPD,CPV,CL,RD,RV,LV0,G,ROWL,EPS,EPSI,CPVMCL PARAMETER (CPD=1005.7, CPV=1870.0, CL=2500.0, RD=287.04) PARAMETER (RV=461.5, LV0=2.501E6, G=9.8, ROWL=1000.0 ) PARAMETER (EPS=RD/RV, EPSI=1./EPS, CPVMCL=CL-CPV) C-------------------------------------------------------------------- C Microphysical parameters: C (here, we use the same values as in the convection scheme) C REAL TLCRIT,ELCRIT,EPMAX PARAMETER (TLCRIT=-55.0, ELCRIT=0.0011, EPMAX=0.999) C-------------------------------------------------------------------- C Initialize output arrays C PRADJ=0.0 C DO I = 1, ND FRADJ(I) = 0.0 FTADJ(I) = 0.0 TNEWLS(I) = T(I) RNEWLS(I) = R(I) QLSP(I) = 0.0 EPLS(I) = 0.0 BIDON(I) = 0.0 ! used for all-or-nothing only SIGSUB(I) = 0.0 ENDDO C------------------------------------------------------------------- C If the cloud parameterization has a "convective" type of closure, C then compute precipitation efficiencies associated with C large-scale precipitation (independent on the presence of C subgrid-scale variability), and then compute the cloud C fraction and in-cloud water content associated with the C combination of subgrid+large-scale condensation. C C DO 9999 I = ND, 1, -1 C C Calculate large-scale condensation and precipitation: C IF(R(I).GT.RS(I))THEN C C Precipitation efficiencies: C EP=0.0 TCA=T(I)-273.15 IF(TCA.GE.0.0)THEN ELACRIT=ELCRIT ELSE ELACRIT=ELCRIT*(1.0-TCA/TLCRIT) END IF ELACRIT=MAX(ELACRIT,0.0) EP= EPMAX * (1.0-ELACRIT/MAX(R(I)-RS(I),1.0E-8)) EP=MAX(EP,0.0) EP=MIN(EP,EPMAX) EPLS(I) = EP ! large-scale precipitation efficiency QLSP(I) = EP*(R(I)-RS(I)) ! precipitated cloud water C C Adjust temperature and humidity profiles: C TCA=T(I)-273.15 ALV=LV0-CPVMCL*TCA CPN=CPD*(1.-R(I))+CPV*R(I) TNEW=(ALV*(EP*R(I)+RS(I)*(ALV/(RV*T(I))-EP)) : +CPN*T(I))/ 1 (CPN+ALV*ALV*RS(I)/(RV*T(I)*T(I))) RNEW=RS(I)*(1.+(TNEW-T(I))*ALV/(RV*T(I)*T(I))) : + (1.-EP)*(R(I)-RS(I)) TNEWLS(I) = TNEW RNEWLS(I) = RNEW FRADJ(I)=FRADJ(I)+(RNEW-R(I))/DT FTADJ(I)=FTADJ(I)-ALV*(RNEW-R(I))/DT/CPN PRADJ=PRADJ-100.0*(PH(I)-PH(I+1))*(RNEW-R(I))/DT : *1000.0*3600.0*24.0/(ROWL*G) END IF ! R>RS 9999 CONTINUE C------------------------------------------------------------- C Cloud parameterization: C c GNO PDF and convect-closure: CALL CLOUDS_GNO_40(ND,R,RS,QSUBGRID,CLDF,CLDQ) C-------------------------------------------------------------------- C Remove large-scale precipitation from the cloud water content: C DO I = 1, ND CLDQ(I) = CLDQ(I) - QLSP(I) ENDDO RETURN END C SUBROUTINE CLOUDS_GNO_40(ND,R,RS,QSUB,CLDF,CLDQ) IMPLICIT NONE C C=================================================================================== C C CLOUDS_GNO version 4.0 C C Purpose: C -------- C C Parameterization of the cloudiness (cloud amount, cloud water content) C associated with cumulus convection. C C Principle: C ---------- C C This cloud parameterization predicts the cloudiness that is associated with C the presence of condensation within a large-scale domain: this condensation C may be produced at the subgrid-scale by cumulus convection and at the C large-scale by super-saturation. C C IMPORTANT: in the present version of the scheme, the only source of subgrid-scale C condensation that is considered is cumulus convection (condensation associated C with boundary layer turbulence, for instance, is not considered). C C The cloud fraction and the in-cloud water content are predicted by a C statistical approach. The subgrid-scale variability of total water C (vapor + condensed) within the gridbox is described by a generalized C log-normal Probability Distribution Function (PDF) whose mean, variance C and skewness coefficient are predicted. The predictors are: C 1) the local concentration of condensed water that is produced at C the subgrid-scale by convection (output of the convection scheme) C 2) the saturation deficit or excess of the environment C 3) the domain-averaged mixing ratio of total water C Note that we impose the distribution of total water to be bounded by zero. C On the other hand, no upper bound of the distribution is considered in this C version of the scheme. C C If no subgrid-scale condensation occurs within the domain, the scheme C becomes equivalent to an "all-or-nothing" large-scale saturation scheme. C C Inputs: C ------- C C ND----------: Number of vertical levels C R--------ND-: Domain-averaged mixing ratio of total water C RS-------ND-: Mean saturation humidity mixing ratio within the gridbox C QSUB-----ND-: Mixing ratio of condensed water within clouds associated C with SUBGRID-SCALE condensation processes (here, it is C predicted by the convection scheme) C Outputs: C -------- C C CLDF-----ND-: cloud fractional area (0-1) C CLDQ-----ND-: in-cloud mixing ratio of condensed water (kg/kg) C C CALL command: C ------------- C C CALL CLOUDS_GNO(ND,R,RS,QSUBGRID,CLDF,CLDQ) C C Reference: C ---------- C C Bony, S and K A Emanuel, 2001: A parameterization of the cloudiness C associated with cumulus convection; Evaluation using TOGA COARE data. C J. Atmos. Sci., accepted. C C Written by: C ----------- C C Sandrine Bony (MIT & LMD/CNRS; bony@wind.mit.edu) - July 2000 C C Difference with version 1.0: C numerical method of resolution of equation 9 C version 1.0: use a Gaussian PDF when erf(v)->1 C version 2.0: use an asymptotic expression of erf(v) instead of a Gaussian PDF C===================================================================================== c c -- input/output arguments of the subroutine: INTEGER ND REAL R(ND), RS(ND), QSUB(ND), CLDF(ND), CLDQ(ND) REAL alpha, lambda, kew, skew, sigs c -- lower bound of the PDF of total water: c REAL PB PARAMETER ( PB = 0.0 ) c c -- parameters controlling the iteration: c -- nmax : maximum nb of iterations (hopefully never reached!) c -- epsilon : accuracy of the numerical resolution (here 2.0%) c -- vmax : v-value above which we use an asymptotic expression for ERF(v) INTEGER nmax, niter PARAMETER ( nmax = 10) REAL epsilon, vmax0, vmax PARAMETER ( epsilon = 0.02, vmax0 = 2.0 ) c -- gardes-fou: REAL min_mu, min_Q PARAMETER ( min_mu = 1.e-12, min_Q=1.e-12 ) c -- misc: INTEGER K, n, m REAL*8 mu, qsat, delta, beta REAL*8 xx, aux, coeff, block, dist, fprime, det REAL*8 pi, u, v, erfu, erfv, xx1, xx2, erfg LOGICAL lconv c ---------------------------------------------------------------------------------- c ---------------------------------------------------------------------------------- pi = ACOS(-1.) c c -- loop over vertical levels : c DO 500 K = 1, ND c mu = R(K) mu = MAX(mu,min_mu) qsat = RS(K) qsat = MAX(qsat,min_mu) delta = log(mu/qsat) IF ( QSUB(K) .lt. min_Q ) THEN C=========================================================================== C If no condensation is produced at the subgrid-scale: C C -> the scheme becomes equivalent to a "large-scale condensation scheme" C ie: cldf = H(mu-qsat) and cldq = (mu-qsat)*H(mu-qsat) C where H is the Heaviside function. C (in the absence of subgrid-scale condensation, the generalized C log-normal PDF becomes equivalent to a gaussian PDF of variance C zero, i.e. it becomes equivalent to a Dirac function and the C cumulative distribution function becomes an Heaviside function). C=========================================================================== CLDQ(K) = MAX( 0.0, mu-qsat ) CLDF(K) = CLDQ(K) / MAX( CLDQ(K), min_mu ) lambda = mu alpha = 0.0 kew = 0.0 skew = 0.0 sigs = 0.0 ELSE C=========================================================================== C Some condensation is produced at the subgrid-scale: C (presence of subgrid-scale variability): C C Use the (iterative) numerical method of Newton to determine the parameters C that characterize the PDF of total water. C C Remark 1: the accuracy of the resolution is controlled by "epsilon" C Remark 2: in GCMs, this numerical method may be too much CPU-time consuming. C In that case, it may be more appropriate to substitute it by a tabulation C of equations 9 and 11 (see the Bony-Emanuel article cited in introduction). C=========================================================================== lconv = .FALSE. ! flag for numerical convergence niter = 0 c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c PDF = generalized log-normal distribution (GNO): c (k<0 if a lower bound is considered for the PDF of total water) c c -> determine x (the parameter k of the GNO PDF) c such that the contribution of subgrid-scale processes to the c in-cloud water content is equal to QSUB(K) c c NB: the "error function" is called ERF or DERF (in double precision) c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ beta = QSUB(K)/mu + EXP( -MIN(0.0,delta) ) vmax = vmax0 IF ( .NOT. lconv ) then ! new c -- roots of equation v > vmax: det = delta + vmax**2. if (det.LE.0.0) vmax = vmax0 + 1.0 det = delta + vmax**2. if (det.LE.0.) then xx = -0.0001 else xx1 = -SQRT(2.)*vmax*(1.0-SQRT(1.0+delta/(vmax**2.))) xx2 = -SQRT(2.)*vmax*(1.0+SQRT(1.0+delta/(vmax**2.))) xx = 1.01 * xx1 if ( xx1 .GE. 0.0 ) xx = 0.5*xx2 endif if (delta.LT.0.) : xx = -0.5*SQRT(log(2.)) ! test comme avant DO n = 1, nmax ! iteration loop u = delta/(xx*sqrt(2.)) + xx/(2.*sqrt(2.)) v = delta/(xx*sqrt(2.)) - xx/(2.*sqrt(2.)) IF ( v .GT. vmax ) THEN IF ( ABS(u) .GT. vmax : .AND. delta .LT. 0. ) THEN c -- use asymptotic expression of erf for u and v large: c ( -> analytic solution for xx ) aux = 2.0*delta*(1.-beta*EXP(delta)) : /(1.+beta*EXP(delta)) xx = -SQRT(aux) block = EXP(-v*v) / v / sqrt(pi) dist = 0.0 fprime = 1.0 ELSE c -- erfv -> 1.0, use an asymptotic expression of erfv for v large: erfu = ERFG(u) aux = sqrt(pi) * (1.0-erfu) * EXP(v*v) coeff = 1.0 - 1./2./(v**2.) + 3./4./(v**4.) block = coeff * EXP(-v*v) / v / sqrt(pi) dist = v * aux / coeff - beta fprime = 2.0 / xx * (v**2.) : * ( coeff*EXP(-delta) - u * aux ) : / coeff / coeff ENDIF ! ABS(u) ELSE c -- general case: erfu = ERFG(u) erfv = ERFG(v) block = 1.0-erfv dist = (1.0 - erfu) / (1.0 - erfv) - beta fprime = 2. /sqrt(pi) /xx /(1.0-erfv)**2. : * ( (1.0-erfv)*v*EXP(-u*u) : - (1.0-erfu)*u*EXP(-v*v) ) ENDIF ! x c -- numerical convergence reached? if ( ABS(dist/beta) .LT. epsilon ) then ! convergence lconv = .TRUE. ! numerical convergence reached with GNO PDF ccc write(*,*) 'CV NEWTON GNO FOR K, N = ',K,n c parameters that characterize the PDF: kew = xx aux = EXP(kew*kew) sigs = mu * SQRT( aux-1.0 ) sigs = MAX(sigs,0.0) lambda = mu / SQRT (aux) alpha = - lambda * kew skew = ( 2.0 + aux**3. - 3*aux ) * (mu/sigs)**3. c deduce the cloud fraction and the in-cloud water content: c -------------------------------------------------------------- CLDF(K) = 0.5 * block CLDQ(K) = QSUB(K) + MAX(mu-qsat,0.0) GOTO 100 else xx = xx - dist/fprime endif ENDDO ! n 100 continue if (.NOT. lconv) then write(*,*) 'NO CV in CLOUDS_GNO: K,mu,qsub,qsat,n,niter: ' : ,K,mu*1000.,QSUB(K)*1000.,qsat*1000.,n,niter write(*,*) 'gamma,beta,error = ',delta,beta,ABS(dist/beta) c all-or-nothing scheme in that (exceptional) case, may be improved later on: CLDQ(K) = MAX( 0.0, mu-qsat ) CLDF(K) = CLDQ(K) / MAX( CLDQ(K), min_mu ) endif ENDIF ! lconv ENDIF ! qsub 500 CONTINUE ! K RETURN END C function erfc(x) parameter( & pv= 9.20888710e+00, & ph= 5.07254732e+00, & p0= 3.86642217e-01, & p1= 1.52430177e-01, & p2= 2.38149125e-02, & p3= 1.30227291e-03, & q0= 1.16381965e-01, & q1= 1.04753802e+00, & q2= 2.92132156e+00, & q3= 6.02608434e+00) y=x*x y=exp(-y)*x*(p3/(y+q3)+p2/(y+q2) & +p1/(y+q1)+p0/(y+q0)) if(x.lt.ph) y=y+2/(exp(pv*x)+1) erfc=y end ! ! error function in double precision ! function erfg(x) implicit real*8 (a - h, o - z) dimension a(0 : 64), b(0 : 64) data (a(i), i = 0, 12) / & 0.00000000005958930743d0, -0.00000000113739022964d0, & 0.00000001466005199839d0, -0.00000016350354461960d0, & 0.00000164610044809620d0, -0.00001492559551950604d0, & 0.00012055331122299265d0, -0.00085483269811296660d0, & 0.00522397762482322257d0, -0.02686617064507733420d0, & 0.11283791670954881569d0, -0.37612638903183748117d0, & 1.12837916709551257377d0 / data (a(i), i = 13, 25) / & 0.00000000002372510631d0, -0.00000000045493253732d0, & 0.00000000590362766598d0, -0.00000006642090827576d0, & 0.00000067595634268133d0, -0.00000621188515924000d0, & 0.00005103883009709690d0, -0.00037015410692956173d0, & 0.00233307631218880978d0, -0.01254988477182192210d0, & 0.05657061146827041994d0, -0.21379664776456006580d0, & 0.84270079294971486929d0 / data (a(i), i = 26, 38) / & 0.00000000000949905026d0, -0.00000000018310229805d0, & 0.00000000239463074000d0, -0.00000002721444369609d0, & 0.00000028045522331686d0, -0.00000261830022482897d0, & 0.00002195455056768781d0, -0.00016358986921372656d0, & 0.00107052153564110318d0, -0.00608284718113590151d0, & 0.02986978465246258244d0, -0.13055593046562267625d0, & 0.67493323603965504676d0 / data (a(i), i = 39, 51) / & 0.00000000000382722073d0, -0.00000000007421598602d0, & 0.00000000097930574080d0, -0.00000001126008898854d0, & 0.00000011775134830784d0, -0.00000111992758382650d0, & 0.00000962023443095201d0, -0.00007404402135070773d0, & 0.00050689993654144881d0, -0.00307553051439272889d0, & 0.01668977892553165586d0, -0.08548534594781312114d0, & 0.56909076642393639985d0 / data (a(i), i = 52, 64) / & 0.00000000000155296588d0, -0.00000000003032205868d0, & 0.00000000040424830707d0, -0.00000000471135111493d0, & 0.00000005011915876293d0, -0.00000048722516178974d0, & 0.00000430683284629395d0, -0.00003445026145385764d0, & 0.00024879276133931664d0, -0.00162940941748079288d0, & 0.00988786373932350462d0, -0.05962426839442303805d0, & 0.49766113250947636708d0 / data (b(i), i = 0, 12) / & -0.00000000029734388465d0, 0.00000000269776334046d0, & -0.00000000640788827665d0, -0.00000001667820132100d0, & -0.00000021854388148686d0, 0.00000266246030457984d0, & 0.00001612722157047886d0, -0.00025616361025506629d0, & 0.00015380842432375365d0, 0.00815533022524927908d0, & -0.01402283663896319337d0, -0.19746892495383021487d0, & 0.71511720328842845913d0 / data (b(i), i = 13, 25) / & -0.00000000001951073787d0, -0.00000000032302692214d0, & 0.00000000522461866919d0, 0.00000000342940918551d0, & -0.00000035772874310272d0, 0.00000019999935792654d0, & 0.00002687044575042908d0, -0.00011843240273775776d0, & -0.00080991728956032271d0, 0.00661062970502241174d0, & 0.00909530922354827295d0, -0.20160072778491013140d0, & 0.51169696718727644908d0 / data (b(i), i = 26, 38) / & 0.00000000003147682272d0, -0.00000000048465972408d0, & 0.00000000063675740242d0, 0.00000003377623323271d0, & -0.00000015451139637086d0, -0.00000203340624738438d0, & 0.00001947204525295057d0, 0.00002854147231653228d0, & -0.00101565063152200272d0, 0.00271187003520095655d0, & 0.02328095035422810727d0, -0.16725021123116877197d0, & 0.32490054966649436974d0 / data (b(i), i = 39, 51) / & 0.00000000002319363370d0, -0.00000000006303206648d0, & -0.00000000264888267434d0, 0.00000002050708040581d0, & 0.00000011371857327578d0, -0.00000211211337219663d0, & 0.00000368797328322935d0, 0.00009823686253424796d0, & -0.00065860243990455368d0, -0.00075285814895230877d0, & 0.02585434424202960464d0, -0.11637092784486193258d0, & 0.18267336775296612024d0 / data (b(i), i = 52, 64) / & -0.00000000000367789363d0, 0.00000000020876046746d0, & -0.00000000193319027226d0, -0.00000000435953392472d0, & 0.00000018006992266137d0, -0.00000078441223763969d0, & -0.00000675407647949153d0, 0.00008428418334440096d0, & -0.00017604388937031815d0, -0.00239729611435071610d0, & 0.02064129023876022970d0, -0.06905562880005864105d0, & 0.09084526782065478489d0 / w = abs(x) if (w .lt. 2.2d0) then t = w * w k = int(t) t = t - k k = k * 13 y = ((((((((((((a(k) * t + a(k + 1)) * t + & a(k + 2)) * t + a(k + 3)) * t + a(k + 4)) * t + & a(k + 5)) * t + a(k + 6)) * t + a(k + 7)) * t + & a(k + 8)) * t + a(k + 9)) * t + a(k + 10)) * t + & a(k + 11)) * t + a(k + 12)) * w else if (w .lt. 6.9d0) then k = int(w) t = w - k k = 13 * (k - 2) y = (((((((((((b(k) * t + b(k + 1)) * t + & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & b(k + 11)) * t + b(k + 12) y = y * y y = y * y y = y * y y = 1 - y * y else y = 1 end if if (x .lt. 0) y = -y erfg = y end ! C SUBROUTINE OPTICAL(ND,T,PH,P,CLDF,CLDQ : ,CLDEMI,CLDTAU,CLDFICE,CLDT,CLDWP) implicit none C============================================================================ C Purpose: C ----------- C Compute cloud optical properties (emissivity, optical thickness) C C CALL command: C ------------- C CALL OPTICAL(NA,T2,PH,P,CLDF,CLDQ,CLDEMI,CLDTAU,CLDFICE,CLDT,CLDWP) C C Inputs: C ----------- C ND-----------: Number of vertical levels C T--------ND--: Grid-box average temperature (K) C PH-------ND+1: Pressure at interface levels (Pa) C P--------ND--: Pressure at mid-levels (Pa) C CLDF-----ND--: Cloud fraction at each vertical level (0-1) C CLDQ-----ND--: In-cloud condensate mixing ratio (kg/kg) C C Outputs: C ----------- C CLDEMI---ND--: Cloud longwave emissivity C CLDTAU---ND--: Cloud optical thickness C CLDFICE--ND--: Ice fraction within the grid-box C CLDT------1--: Total cloud cover C CLDWP-----1--: Gridbox averaged cloud water path (kg/m2) C C Author & date: C --------------- C Sandrine Bony (MIT & LMD/CNRS) - June 1999 C C Modification Feb 2000: if newmicro=TRUE, use cloud microphysical properties C suggested by Iacobellis and Somerville (1999) for C rad_liq and rad_ice, and by Ebert and Curry (1992) C for coef_liq and coef_ice. C C Caution: if newmicro=TRUE and if this routine is used in a GCM, information C about the surface (land or ocean) should be used to specify the C droplet concentration (Noc or Nland); for the moment we C assume that this is an ocean surface. C C============================================================================ C integer ND,k,i,n,no real T(ND),CLDF(ND),CLDQ(ND),CLDEMI(ND),CLDTAU(ND) : ,PH(ND+1),P(ND),CLDFICE(ND),CLDT,CLDWP c c -- equivalent cloud droplet radius for liquid/ice clouds: real radius, rad_ice, rad_liq parameter (rad_liq=10.0, rad_ice=30.0) c c -- IR absorption coefficient for liquid/ice clouds (includes c the diffusivity factor): real coef, coef_ice, coef_liq parameter (coef_liq=0.13, coef_ice=0.09) c c -- Thresholds for liquid/ice clouds distinction: real seuil_neb, T_ice, T_o parameter (seuil_neb=0.001) c from LMD ("CTRL" and "ICE-OPT" expts): parameter (T_o=273.15, T_ice=273.15-15.0) c Zender and Kiehl (95): ("PHASE" expt) c parameter (T_o=273.15-10, T_ice=273.15-30.0) c c -- for the new microphysical parameterization: real Noc, Nland ! droplet nb concentration (cm^{-3}) parameter (Noc=150., Nland=600.)! Bower et al. 1994 real k_liq, k_ice0, k_ice, DF, k0,k1,k2,k3 parameter (k_liq=0.0903, k_ice0=0.005) ! units = m2/g parameter (DF=1.66) ! diffusivity factor parameter (k0=60.75, k1=-2.47, k2=-0.11, k3=-0.001) real pi, aux, rel, rei, kabs, tc c -- Misc: real RG, undef, zflwp, zfiwp parameter (RG=9.8,undef=999.999) logical lo, newmicro parameter (newmicro=.FALSE.) !if T: new opt prop for ice clds, refnew c parameter (newmicro=.TRUE.) !if T: new opt prop for ice clds pi = ACOS(-1.) c c Cloud optical thickness and emissivity: c DO k = 1, ND CLDF(k) = MAX(CLDF(k), seuil_neb) c liquid/ice fraction: CLDFICE(k) = 1.0 - (T(k)-T_ice) / (T_o-T_ice) CLDFICE(k) = MIN(MAX(CLDFICE(k),0.0),1.0) c liquid/ice cloud water paths: zflwp = 1000.*(1.-CLDFICE(k))*CLDQ(k) : *(PH(k)-PH(k+1))/RG zfiwp = 1000.*CLDFICE(k)*CLDQ(k) : *(PH(k)-PH(k+1))/RG c optical properties: IF (newmicro) THEN c*********************************************************** c NB: "no" means that we modify only the optical properties c of ice clouds, not those of liquid clouds c*********************************************************** c -- parameterization of the effective cloud droplet radius (microns): c++++ for liquid water clouds: rel = rad_liq ! no c++++ for ice clouds: as a function of the ambiant temperature tc = T(k)-273.15 c ...... formula used by Iacobellis and Somerville (2000): c ...... (with an asymptotical value of 3.5 microns at T<-81.4 C c ...... added to be consistent with observations of Heymsfield et al. 1986): rei = 0.71*tc + 61.29 if (tc.le.-81.4) rei = 3.5 ! only micronew, phasenew expts c cloud optical thickness: c for liquid clouds, LMD-like, for ice clouds, Ebert & Curry (1992): c ------------------------------------------------------------------ if (zflwp.eq.0.) rel = 1. ! no influence if (zfiwp.eq.0. .or. rei.le.0.) rei = 1. ! no influence CLDTAU(k) = 3.0/2.0 * ( zflwp/rel ) . + zfiwp * (3.448e-03 + 2.431/rei) c cloud infrared emissivity: c --------------------------- c the broadband infrared absorption coefficient is parameterized as a c function of the effective cld droplet radius: c ... Ebert and Curry (1992) formula: c as used by Kiehl & Zender (1995): k_ice = k_ice0 + 1.0/rei ! iceopt c cloud emissivity: CLDEMI(k) = 1.0 . - EXP( - coef_liq*zflwp - DF*k_ice*zfiwp ) ELSE c -- version LMD: c --------------- CLDTAU(k) = 3.0/2.0 . * ( zflwp/rad_liq + zfiwp/rad_ice ) CLDEMI(k) = 1.0 . - EXP( - coef_liq*zflwp - coef_ice*zfiwp ) ENDIF ! newmicro lo = (CLDF(k) .LE. seuil_neb) IF (lo) CLDF(k) = 0.0 IF (lo) CLDTAU(k) = 0.0 IF (lo) CLDEMI(k) = 0.0 lo = (CLDQ(k) .EQ. undef) IF (lo) write(*,*) 'PB: CLDQ EQ UNDEF ' IF (lo) CLDTAU(k) = undef IF (lo) CLDEMI(k) = undef ENDDO C C Cloud liquid path and total cloudiness: C CLDT = 1. CLDWP = 0. no = 0 DO k = ND, 1, -1 if (CLDQ(k) .NE. undef .and. CLDF(k).ne.undef) then CLDWP = CLDWP . + CLDQ(k)*CLDF(k)*(PH(k)-PH(k+1))/RG CLDT = CLDT * (1.0-CLDF(k)) ! random overlap else no = 1 write(*,*) 'PB: k, CLDF, CLDQ: ',CLDF(k),CLDQ(k) endif ENDDO CLDT = 1.0 - CLDT ! random overlap only C if (no.eq.1) then CLDT = undef CLDWP = undef endif RETURN END c====================================================================== SUBROUTINE orbite(xjour,longi,dist) IMPLICIT none c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) (adapte du GCM du LMD) date: 19930818 c Objet: pour un jour donne, calculer la longitude vraie de la terre c (par rapport au point vernal-21 mars) dans son orbite solaire c calculer aussi la distance terre-soleil (unite astronomique) c====================================================================== c Arguments: c xjour--INPUT--R- jour de l'annee a compter du 1er janvier c longi--OUTPUT-R- longitude vraie en degres par rapport au point c vernal (21 mars) en degres c dist---OUTPUT-R- distance terre-soleil (par rapport a la moyenne) REAL xjour, longi, dist c====================================================================== include "YOMCST.h" C C -- Variables dynamiques locales REAL pir,xl,xllp,xee,xse,xlam,dlamm,anm,ranm,anv,ranv C C -- sb: c call suphec to fill YOMSCT: call suphec C sb -- pir = 4.0*ATAN(1.0) / 180.0 xl=R_peri+180.0 xllp=xl*pir xee=R_ecc*R_ecc xse=SQRT(1.0-xee) xlam = (R_ecc/2.0+R_ecc*xee/8.0)*(1.0+xse)*SIN(xllp) . - xee/4.0*(0.5+xse)*SIN(2.0*xllp) . + R_ecc*xee/8.0*(1.0/3.0+xse)*SIN(3.0*xllp) xlam=2.0*xlam/pir dlamm=xlam+(xjour-81.0) anm=dlamm-xl ranm=anm*pir xee=xee*R_ecc ranv=ranm+(2.0*R_ecc-xee/4.0)*SIN(ranm) . +5.0/4.0*R_ecc*R_ecc*SIN(2.0*ranm) . +13.0/12.0*xee*SIN(3.0*ranm) c anv=ranv/pir longi=anv+xl C dist = (1-R_ecc*R_ecc) . /(1+R_ecc*COS(pir*(longi-(R_peri+180.0)))) RETURN END c====================================================================== SUBROUTINE angle(longi, lati, frac, muzero) IMPLICIT none c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818 c Objet: Calculer la duree d'ensoleillement pour un jour et la hauteur c du soleil (cosinus de l'angle zinithal) moyenne sur la journee c====================================================================== c Arguments: c longi----INPUT-R- la longitude vraie de la terre dans son plan c solaire a partir de l'equinoxe de printemps (degre) c lati-----INPUT-R- la latitude d'un point sur la terre (degre) c frac-----OUTPUT-R la duree d'ensoleillement dans la journee divisee c par 24 heures (unite en fraction de 0 a 1) c muzero---OUTPUT-R la moyenne du cosinus de l'angle zinithal sur c la journee (0 a 1) c====================================================================== include "dimensions.h" include "dimphy.h" REAL longi REAL lati(klon), frac(klon), muzero(klon) include "YOMCST.h" REAL lat, omega, lon_sun, lat_sun REAL pi_local, incl INTEGER i c pi_local = 4.0 * ATAN(1.0) incl=R_incl * pi_local / 180. c lon_sun = longi * pi_local / 180.0 lat_sun = ASIN (sin(lon_sun)*SIN(incl) ) c DO i = 1, klon lat = lati(i) * pi_local / 180.0 lat = lati(i) * pi_local / 180.0 c IF ( lat .GE. (pi_local/2.+lat_sun) . .OR. lat.LE.(-pi_local/2.+lat_sun)) THEN omega = 0.0 ! nuit polaire ELSE IF ( lat.GE.(pi_local/2.-lat_sun) . .OR. lat.LE.(-pi_local/2.-lat_sun)) THEN omega = pi_local ! journee polaire ELSE omega = -TAN(lat)*TAN(lat_sun) omega = ACOS (omega) ENDIF c frac(i) = omega / pi_local c IF (omega .GT. 0.0) THEN muzero(i) = SIN(lat)*SIN(lat_sun) . + COS(lat)*COS(lat_sun)*SIN(omega) / omega ELSE muzero(i) = 0.0 ENDIF ENDDO c RETURN END c==================================================================== SUBROUTINE zenang(longi,gmtime,pdtrad,lat,long, s pmu0,frac) IMPLICIT none c============================================================= c Auteur : O. Boucher (LMD/CNRS) c d'apres les routines zenith et angle de Z.X. Li c Objet : calculer les valeurs moyennes du cos de l'angle zenithal c et l'ensoleillement moyen entre gmtime1 et gmtime2 c connaissant la declinaison, la latitude et la longitude. c Rque : Different de la routine angle en ce sens que zenang c fournit des moyennes de pmu0 et non des valeurs c instantanees, du coup frac prend toutes les valeurs c entre 0 et 1. c Date : premiere version le 13 decembre 1994 c revu pour GCM le 30 septembre 1996 c=============================================================== c longi----INPUT : la longitude vraie de la terre dans son plan c solaire a partir de l'equinoxe de printemps (degre) c gmtime---INPUT : temps universel en fraction de jour c pdtrad---INPUT : pas de temps du rayonnement (secondes) c lat------INPUT : latitude en degres c long-----INPUT : longitude en degres c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad c================================================================ include "dimensions.h" include "dimphy.h" include "YOMCST.h" c================================================================ real longi, gmtime, pdtrad real lat(klon), long(klon), pmu0(klon), frac(klon) c================================================================ integer i real gmtime1, gmtime2 real pi_local, deux_pi_local, incl real omega1, omega2, omega c omega1, omega2 : temps 1 et 2 exprime en radian avec 0 a midi. c omega : heure en radian du coucher de soleil c -omega est donc l'heure en radian de lever