subroutine O3Merra_long_60 c c rearranged for GOME data 24.03.2007, Alex Pogoreltsev c c----------------- nb=36 is fixed for given O3 model c -----------------and kgit has to be set (arbitrary) c--------------------------------------------------------------------- include 'com_main.fc' parameter (klev=23) c------------------------------ parameter (nb=36,kgit=48,igit=64) real o3vmr (klev,nb,igit),blev(klev),o3out2(nb,klev) real rlatO3(nb),dolgO3(igit) common/o3_Berl/blev,o3out2,h_km c -------------------------------------------------------------------- c------ CO3(45,nb) will be used to calculate O3 cooling c c O3ppmv - O3 part per million volume c CO3 - for O3 cooling c SO3 - Integral of O3 number density c VMOINT - integralnaya koncentraciya O3 c VMO - koncentraciya O3 c--------------------------------------------------------------------- EXTERNAL VMOINT REAL Z_km(kgit+1) h_km=7. c-------------------------------- READ O3 data open(4,file='oz_merra_03_0414.dat') Do k=1,klev c-------------------------- zkm_m(k)=-7.0*alog(egb(k)/1000.0) do j=1,nb do i=1,igit read(4,4) blev(k),rlatO3(j),dolgO3(i),o3vmr(k,j,i) 4 format(1x,3f8.2,e13.5) end do end do End do close (4) do j=1,nb print*,o3vmr(5,j,27),blev(5), nb,klev,igit,kgit end do c------------------ O3_vmr averaged over longitude c do k=1,klev c do j=1,nb c o3out2(j,k)=0. c DO 1111 i=1,igit c o3out2(j,k)=o3out2(j,k)+o3vmr(k,j,i)/float(igit) c 1111 CONTINUE c if(o3out2(j,k).le.0.) print*, j,k c end do c end do c print*,o3out2(18,10) c----------------------------------------------------------------- DO 1112 i=1,igit c------------------ version with longitudinal variations of O3_vmr do k=1,klev do j=1,nb o3out2(j,k)=o3vmr(k,j,i) if(o3out2(j,k).le.0.) o3out2(j,k)=1.e-6 end do end do print*,o3out2(18,10) c----------------------------------------------------------------- DZKM=135./(48.-0.5) ZUNTKM=0.5*DZKM Z_km(1)=0. DO 20 j=1,NB DO 10 K=1,KGIT Z_km(K+1)=ZUNTKM+FLOAT(K-1)*DZKM Z_out=Z_km(k+1) OZPPMV(j,K,i)=VMO(Z_out,j) 10 CONTINUE c print*,ozppmv(23,5,27) c------------------------------ O3 ppmv for cooling DO 11 K=1,45 Z_out=FLOAT(K-1)*0.25*h_km CO3(K,j,i)=VMO(Z_out,j) 11 CONTINUE 20 CONTINUE c---------------------------------------------------------------- ZOBEN1=Z_km(11) ZOBEN2=200. EPSMAX=1.E-6 DO 42 j=1,NB DO 41 K=1,KGIT IF(K.LT.10) THEN CALL ADAPT(j,VMOINT,Z_km(K+1),ZOBEN1,EPS,EPSMAX,RES1) CALL ADAPT(j,VMOINT,ZOBEN1 ,ZOBEN2,EPS,EPSMAX,RES2) RES=RES1+RES2 ELSE CALL ADAPT(j,VMOINT,Z_km(K+1),ZOBEN2,EPS,EPSMAX,RES) END IF c CALL ADAPT(PHI(j),j,VMOINT,Z_km(K+1),ZOBEN,EPS,EPSMAX,SW) c-------------- N_O3*1e5(km->cm)*P_0/K_B/T_0 SO3(j,K,i)=RES*1.e5*1.013/1.38066e-16/273. 41 CONTINUE 42 CONTINUE 1112 CONTINUE print*, nb C RETURN END C-------------------------------------------------------------------------- C SUBROUTINE ADAPT (II,FCT,A1,B1,EPS,EPSMAX,XINT) EXTERNAL FCT DIMENSION QI(25),BI(25) NMAX = 25 XFF = 8.192E3 XFAK = 1./(XFF - 1.) C BI(1) = B1 XA = A1 CALL GAUSS (II,FCT,QI(1),XA,B1) C XINT = 0. EPS = 0. J = 1 C 100 IF (J .LT. 1) RETURN QI1 = QI(J) BIE = BI(J) XAP = 0.5*(XA + BIE) CALL GAUSS (II,FCT,QI2,XA,XAP) C CALL GAUSS (II,FCT,QI3,XAP,BIE) C QI0 = QI2 + QI3 DELT = XFAK*ABS(QI0 - QI1) DREL = DELT IF (QI0 .NE. 0.) DREL = DELT/ABS(QI0) IF (DREL .LE. EPSMAX) GO TO 10 C JP1 = J + 1 IF (JP1 .GT. NMAX) GO TO 10 QI(J) = QI3 C J = JP1 BI(J) = XAP QI(J) = QI2 GO TO 100 C 10 XINT = XFAK*(XFF*QI0 - QI1) + XINT EPS = EPS + DELT XA = BIE J = J - 1 GO TO 100 C END C C----------------------------------------------------------- C SUBROUTINE GAUSS (II,FCT,VALUE,AF,BF) C DIMENSION STUETZ(6) , GEW(6) DATA STUETZ /-0.9324695,-0.6612094,-0.2386192,0.2386192 *,0.6612094,0.9324695/ DATA GEW /0.1713245,0.3607616,0.4679139,0.4679139, * 0.3607616,0.1713245/ C HH1 = 0.5*(BF - AF) VALUE = 0. XXM = 0.5*(AF + BF) DO 100 I = 1 , 6 VAR = XXM + HH1*STUETZ(I) SUM = FCT (VAR,II)*GEW(I) 100 VALUE = VALUE + SUM C VALUE = HH1*VALUE RETURN END C C----------------------------------------------------------- FUNCTION VMOINT(ZZ,j) Parameter(klev=23,nb=36) Real blev(klev),o3lev(klev),o3out2(nb,klev), $ B(klev),C(klev),D(klev) common/o3_Berl/blev,o3out2,h_km do k=1,klev o3lev(k)=ALOG(o3out2(j,k)) end do CALL SPL3(klev,blev,o3lev,B,C,D) if(zz.le.blev(klev)) then CALL SEVAL0(klev,zz,blev,o3lev,B,C,D,O3ppmv) O3ppmv=EXP(O3ppmv) if(O3ppmv.le.1.e-6) O3ppmv=1.e-6 else O3ppmv=EXP(o3lev(klev))*exp((blev(klev)-zz)/h_km) end if VMOINT=O3ppmv*EXP(-ZZ/h_km) RETURN END C-------------------------------------------------------------- FUNCTION VMO(ZZ,j) Parameter(klev=23,nb=36) real blev(klev),o3lev(klev),o3out2(nb,klev), $ B(klev),C(klev),D(klev) common/o3_Berl/blev,o3out2,h_km do k=1,klev o3lev(k)=ALOG(o3out2(j,k)) end do CALL SPL3(klev,blev,o3lev,B,C,D) if(zz.le.blev(klev)) then CALL SEVAL0(klev,zz,blev,o3lev,B,C,D,O3ppmv) O3ppmv=EXP(O3ppmv) else O3ppmv=EXP(o3lev(klev))*exp((blev(klev)-zz)/h_km) end if VMO=O3ppmv RETURN END c-------------------------------------------------------------- SUBROUTINE SPL3(N,X,Y,B,C,D) REAL T, X(N), Y(N), B(N), C(N), D(N) NM1 = N - 1 IF( N .LT. 2 ) RETURN IF( N .LT. 3 ) GO TO 50 D(1) = X(2)-X(1) C(2) = (Y(2)-Y(1))/D(1) DO 10 I = 2, NM1 D(I) = X(I+1)-X(I) B(I) = 2.*(D(I-1)+D(I)) C(I+1) = (Y(I+1)-Y(I))/D(I) C(I) = C(I+1)-C(I) 10 CONTINUE B(1) = -D(1) B(N) = -D(N-1) C(1) = 0. C(N) = 0. IF( N .EQ. 3 ) GO TO 15 C(1) = C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1)) C(N) = C(N-1)/(X(N)-X(N-2)) * -C(N-2)/(X(N-1)-X(N-3)) C(1) = C(1)*D(1)**2/(X(4)-X(1)) C(N) =-C(N)*D(N-1)**2/(X(N)-X(N-3)) 15 CONTINUE DO 20 I = 2, N T = D(I-1)/B(I-1) B(I) = B(I)-T*D(I-1) C(I) = C(I)-T*C(I-1) 20 CONTINUE C(N) = C(N)/B(N) DO 30 IB = 1, NM1 I = N - IB C(I) = (C(I)-D(I)*C(I+1))/B(I) 30 CONTINUE B(N) = (Y(N)-Y(NM1))/D(NM1) * + D(NM1)*(C(NM1)+2.*C(N)) DO 40 I = 1, NM1 B(I) = (Y(I+1)-Y(I))/D(I) * -D(I)*(C(I+1)+2.*C(I)) D(I) = (C(I+1)-C(I))/D(I) C(I) = 3.*C(I) 40 CONTINUE C(N) = 3.*C(N) D(N) = D(N-1) RETURN 50 CONTINUE B(1) = (Y(2)-Y(1))/(X(2)-X(1)) C(1) = 0. D(1) = 0. B(2) = B(1) C(2) = 0. D(2) = 0. RETURN END c---------------------------------------------------------- SUBROUTINE SEVAL0(N,U,X,Y,B,C,D,S) REAL U, X(N), Y(N), B(N), C(N), D(N), S DATA I/1/ IF( I .GE. N ) I = 1 IF( U .LT. X(I) ) GO TO 10 IF( U .LE. X(I+1) ) GO TO 30 10 CONTINUE I = 1 J = N + 1 20 CONTINUE K = (I+J) / 2 IF( U .LT. X(K) ) J = K IF( U .GE. X(K) ) I = K IF( J .GT. I+1 ) GO TO 20 30 CONTINUE DX = U - X(I) S = Y(I) + DX*(B(I)+DX*(C(I)+DX*D(I))) RETURN END