#include "dims.h" SUBROUTINE CHAPMN use init_module,only: secs,sin_sundec,cos_sundec use cons_module,only: kmaxp1,t0,rmass,expz,len1,len3,kmax,pi, | dphi,dlamda,re,avo,p0,expzmid,grav,gask implicit none C **** CALCULATE LINE INTEGRAL FOR EACH SPECIES AT EACH GRID C **** POINT IN NNO2, NNO, NNN2 #include "params.h" #include "fcom.h" #include "vscr.h" #include "index.h" #include "buff.h" #include "phys.h" ! ! Local: real :: rlat,coslat,sinlat,factor,sqrtpi,exparg integer :: i,nzk,ntk,n,k,nnk,nps1k,nps2k,npsk,nmsk,kk,nnvk real,parameter :: big=1.e80 real :: t1a(zimxp) real,dimension(zimxp,zkmxp) :: ! for diagnostic plotting | chi_plt, slt_plt, idn_plt C **** T1=SOLAR ZENITH ANGLE,CHI C **** T2=SIN(CHI) C **** T3=COS(CHI) C **** T4=SQRT(SIN(CHI)) sqrtpi = sqrt(pi) ! 1/00: removing nhemi ! RLAT=-.5*pi*FLOAT(1-NHEMI)+(FLOAT(J-1)+.5)*dphi RLAT=-.5*pi+(FLOAT(J-1)+.5)*dphi COSLAT=COS(RLAT) SINLAT=SIN(RLAT) DO 23 I=1,LEN1 T1(I)=FLOAT(I-3) 23 CONTINUE DO 1 I=1,LEN1 C **** T1=LOCAL TIME T1a(I)=AMOD(SECS/3600.+(T1(I)*dlamda+pi)*12./pi,24.) C **** T1=CHI T1(I)=ACOS(sin_sundec*SINLAT+cos_sundec*COSLAT*COS(pi* | (T1a(I)-12.)/12.)) T2(I)=SIN(T1(I)) T3(I)=COS(T1(I)) T4(I)=SQRT(T2(I)) chi_plt(i,:) = t1(i) slt_plt(i,:) = t1a(i) 1 CONTINUE ! ! idn(zimxp) is day-night index (day=1, night=0) for current latitude ! (see fcom.h): idn(:) = 1 do i=1,len1 if (t1(i) > 1.8326) idn(i) = 0 idn_plt(i,:) = idn(i) enddo ! call addfsech('CHI',' ',' ',chi_plt,zimxp,zkmxp,zkmxp,j) ! call addfsech('SLT',' ',' ',slt_plt,zimxp,zkmxp,zkmxp,j) ! call addfsech('IDN',' ',' ',idn_plt,zimxp,zkmxp,zkmxp,j) C **** S9=RP NZK=NJ+NZ DO 2 I=1,LEN3 S9(I,1) = F(I,NZK)+re 2 CONTINUE ! call addfsech('RP',' ',' ',s9,zimxp,zkmxp,zkmxp,j) C **** S8=T (INTERPOLATED TO KMAXP1 LEVELS) C **** LEVEL 1 NTK=NJ+NT+KMAX DO 3 I=1,LEN1 S8(I,1)=F(I,NTK)+T0(1) 3 CONTINUE C **** LEVELS 2 THRU KMAX NTK=NJ+NT-1 DO 4 K=2,KMAX NTK=NTK+1 DO 4 I=1,LEN1 S8(I,K)=.5*(F(I,NTK)+F(I,NTK+1))+T0(K) 4 CONTINUE C **** LEVEL KMAXP1 NTK=NJ+NT+KMAX-1 DO 5 I=1,LEN1 S8(I,KMAXP1)=F(I,NTK)+T0(KMAXP1) 5 CONTINUE ! call addfsech('TNI',' ',' ',s8,zimxp,zkmxp,zkmxp,j) C **** SPECIES DO LOOP DO 6 N=1,3 C **** S7=PSI(N) IF(N.EQ.3)THEN NPS1K=NJ+NPS NPS2K=NPS1K+KMAXP1 DO 7 I=1,LEN3 S7(I,1)=1.-F(I,NPS1K)-F(I,NPS2K) 7 CONTINUE ELSE NPSK=NJ+NPS+(N-1)*KMAXP1 DO 8 I=1,LEN3 S7(I,1)=F(I,NPSK) 8 CONTINUE ENDIF do i=1,len3 if (s7(i,1) < 0.) s7(i,1) = 0. enddo ! if (n==1) then ! call addfsech('CHAP_O2',' ',' ',s7,zimxp,zkmxp,zkmxp,j) ! elseif (n==2) then ! call addfsech('CHAP_O1',' ',' ',s7,zimxp,zkmxp,zkmxp,j) ! else ! call addfsech('CHAP_N2',' ',' ',s7,zimxp,zkmxp,zkmxp,j) ! endif C **** S7=NO. DENSITY INTEGRAL C **** TOP OF MODEL C **** S7(KMAXP1)=N0*P0*PSI(N)*expz*MBAR/(RMASS(N)**2*G) FACTOR=avo*p0*expz(KMAX)*expzmid/(RMASS(N)**2*grav) NMSK=NJ+NMS+KMAX DO 9 I=1,LEN1 S7(I,KMAXP1)=FACTOR*.5*(S7(I,KMAX)+S7(I,KMAXP1))*F(I,NMSK) 9 CONTINUE C **** INTEGRATE DOWN TO LEVEL 1 DO 10 KK=1,KMAX K=KMAXP1-KK FACTOR=avo*p0*expz(K)/(RMASS(N)*grav)*dz DO 10 I=1,LEN1 S7(I,K)=S7(I,K+1)+FACTOR*S7(I,K) 10 CONTINUE C **** COPY VERTICAL COLUMN DENSITIES TO NNVO2, NNVO, NNVN2 NNVK=NNVO2+(N-1)*KMAXP1 DO 16 I=1,LEN3 F(I,NNVK)=S7(I,1) 16 CONTINUE ! if (n==1) then ! call addfsech('VO2',' ',' ',f(1,nnvo2),zimxp,zkmxp,zkmx,j) ! elseif (n==2) then ! call addfsech('VO1',' ',' ',f(1,nnvo),zimxp,zkmxp,zkmx,j) ! elseif (n==3) then ! call addfsech('VN2',' ',' ',f(1,nnvn2),zimxp,zkmxp,zkmx,j) ! endif C **** S6=SQRT(RP/2HP) FACTOR=RMASS(N)*grav/(2.*gask) DO 11 I=1,LEN3 S6(I,1)=SQRT(S9(I,1)*FACTOR/S8(I,1)) 11 CONTINUE C **** S5=YP DO 12 K=1,KMAXP1 DO 12 I=1,LEN1 S5(I,K)=S6(I,K)*ABS(T3(I)) 12 CONTINUE C **** S5=IP DO 13 I=1,LEN3 if (s5(i,1) >= 8.) then s5(i,1) = S7(I,1)*sqrtpi*S6(I,1)*(0.56498823/(0.06651874+ | S5(I,1))) else s5(i,1) = S7(I,1)*sqrtpi*S6(I,1)*((1.0606963+0.5564383* | S5(I,1))/((S5(I,1)+1.7245609)*S5(I,1)+1.0619896)) endif 13 CONTINUE C **** S4=2.*IG FACTOR=grav*RMASS(N)/gask DO K=1,KMAXP1 DO I=1,LEN1 ! ! FP overflow here on prospect (probably at exp(x) where x is too big) ! 3/10/00: problem at j=13. Compile with -fpe2 to avoid abort. ! 5/22/01: use conditional on exparg and idn, as in tgcm24. ! S4(I,K)=2.*S7(I,K)*EXP(S9(I,K)*(1.-T2(I))*FACTOR/S8(I,K))* ! | sqrtpi*T4(I)*S6(I,K) ! exparg = S9(I,K)*(1.-T2(I))*FACTOR/S8(I,K) if (idn(i)==1.and.exparg < 650.) then ! daytime S4(I,K)=2.*S7(I,K)*EXP(exparg)*sqrtpi*T4(I)*S6(I,K) else s4(i,k) = big endif enddo enddo C **** SN=LINE INTEGRAL (0.0 IF OBSCURED BY EARTH) NNK=NNO2+(N-1)*KMAXP1-1 DO K=1,KMAXP1 NNK=NNK+1 DO I=1,LEN1 if (t3(i) >= 0.) then f(i,nnk) = 1.e80 else f(i,nnk) = s4(i,k)-s5(i,k) endif if (S9(I,K)*T2(I)-re < 0.) f(i,nnk) = 1.e80 if (t3(i) >= 0.) f(i,nnk) = s5(i,k) enddo ! if (any(f(:,nnk)==big)) then ! write(6,"('chapmn: lat=',i2,' n=',i2,' k=',i2, ! | ' f(:,nnk)=',/,(6e12.4))") j,n,k,f(:,nnk) ! endif enddo ! if (n==1) then ! call addfsech('SCO2',' ',' ',f(1,nno2),zimxp,zkmxp,zkmx,j) ! elseif (n==2) then ! call addfsech('SCO1',' ',' ',f(1,nno),zimxp,zkmxp,zkmx,j) ! elseif (n==3) then ! call addfsech('SCN2',' ',' ',f(1,nnn2),zimxp,zkmxp,zkmx,j) ! endif 6 CONTINUE RETURN END C