! subroutine comp(tn,o2,o2_nm,o1,o1_nm,un,vn,w,hdo2,hdo1, | o2_upd,o2nm_upd,o1_upd,o1nm_upd, | lev0,lev1,lon0,lon1,lat0,lat1) ! ! This software is part of the NCAR TIE-GCM. Use is governed by the ! Open Source Academic Research License Agreement contained in the file ! tiegcmlicense.txt. ! ! Advance major species O2 and O. ! use params_module,only: dz,nlonp4,nlat,dlat use init_module,only: glat,iday use cons_module,only: pi,hor,dtr,rmassinv_o2,rmassinv_o1, | rmassinv_n2,rmass_o2,rmass_o1,expz,expzmid,expzmid_inv, | difk,dtx2inv,dtsmooth,dtsmooth_div2,difhor use chemrates_module,only: fs ! from sub comp_o2o use lbc,only: b,fb use fields_module,only: tlbc use addfld_module,only: addfld use diags_module,only: mkdiag_O_N2 implicit none ! ! Args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1 real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: | tn, ! neutral temperature | o2, ! O2 (mmr) at current timestep | o1, ! O (mmr) at current timestep | o2_nm,! O2 (mmr) at time n-1 | o1_nm,! O (mmr) at time n-1 | un, ! zonal wind velocity at current timestep | vn, ! meridional wind velocity at current timestep | w, ! vertical velocity at current timestep | hdo2, ! O2 horizontal diffusion (hdif3) | hdo1 ! O horizontal diffusion (hdif3) real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2), | intent(out):: | o2_upd ,o1_upd, ! output: O2,O updated for next timestep | o2nm_upd ,o1nm_upd ! output: O2,O updated for previous timestep ! ! Local: integer :: k,kk,i,lat,isp,km,kp,ktmp,m,lonbeg,lonend integer :: nk,nkm1,i0,i1,k0,k1 integer,parameter :: io2=1,io1=2 ! indices to O2, O, respectively real,dimension(lon0:lon1,lev0:lev1,2,2) :: gama real,dimension(lon0:lon1,lev0:lev1,2) :: zz real,dimension(lon0:lon1,lev0:lev1) :: embar real,dimension(lon0:lon1,2,2,2) :: ak real,dimension(lon0:lon1,2,2) :: ep,pk,qk,rk,wkm1,wkm2 real,dimension(lon0:lon1,2) :: fk,wkv1,wkv2,ps0 real,dimension(lon0:lon1) :: wks1,wks2,wks3,wks4, | embar0,dfactor real :: rlat real :: ak0(2,2),phi(2,3),delta(2,2),tau,t00,small real,dimension(lev0:lev1,lon0:lon1,lat0:lat1) :: | o2nm_smooth, o1nm_smooth, ! smoothed at time n-1 | o2_advec , o1_advec ! horizontal advection real,dimension(lev0:lev1,lon0:lon1,2) :: upd ! For diagnostics: real,dimension(lev0:lev1,lon0:lon1,2) :: zz_ki real,dimension(lev0:lev1,lon0:lon1,2,2) :: gama_ki ! phi(:,1)=(/0. ,0.673/) phi(:,2)=(/1.35,0. /) phi(:,3)=(/1.11,0.769/) tau=1.86e+3 delta(:,1)=(/1.,0./) delta(:,2)=(/0.,1./) t00=273. small = 1.e-6 i0=lon0 ; i1=lon1 ; k0=lev0 ; k1=lev1 ! nk = lev1-lev0+1 nkm1 = nk-1 ! ! Calculate and save horizontal advection in o2_advec, o1_advec: ! do lat=lat0,lat1 call advecl(o2,o1,un,vn,o2_advec,o1_advec, | lev0,lev1,lon0,lon1,lat0,lat1,lat) enddo ! lat=lat0,lat1 ! ! Save smoothed o2,o at time n-1: ! call smooth(o2_nm,o2nm_smooth,lev0,lev1,lon0,lon1,lat0,lat1,0) call smooth(o1_nm,o1nm_smooth,lev0,lev1,lon0,lon1,lat0,lat1,0) ! do lat=lat0,lat1 ! call addfld('O2SMOOTH','O2SMOOTH',' ',o2nm_smooth(:,i0:i1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! call addfld('O1SMOOTH','O1SMOOTH',' ',o1nm_smooth(:,i0:i1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! enddo ! lat=lat0,lat1 ! ! Begin latitude scan: do lat=lat0,lat1 ! do i=lon0,lon1 dfactor(i) = 1. enddo ! ! difhor flag is a parameter flag set to 1 in cons.F ! dfactor = ! .5*(1.+SIN(PI*(ABS(RLATM)-PI/6.)/(PI/3.))) FOR ABS(RLATM).LT.PI/3. ! dfactor = 1. FOR ABS(RLATM).GE.PI/3 ! (dfactor was in sub dfact in earlier versions) ! if (difhor > 0) then rlat = (glat(1)+(lat-1)*dlat)*dtr if (abs(rlat)-pi/4.5 >= 0.) then dfactor(:) = hor(lat)+1. else dfactor(:) = hor(lat)+.5*(1.+sin(pi*(abs(rlat)-pi/9.)/ | (pi/4.5))) endif else dfactor(:) = 1. endif ! write(6,"('comp: lat=',i2,' dfactor=',/,(6e12.4))") ! | lat,dfactor ! ! Embar: do i=lon0,lon1 do k=lev0,lev1 embar(i,k) = 1./(o2(k,i,lat)*rmassinv_o2 + | o1(k,i,lat)*rmassinv_o1 + | (1.-o2(k,i,lat)-o1(k,i,lat))*rmassinv_n2) enddo ! k=lev0,lev1 enddo ! i=lon0,lon1 ! real,dimension(lon0:lon1,lev0:lev1) :: embar ! call addfld('EMBAR',' ',' ',embar(i0:i1,:), ! | 'lon',i0,i1,'lev',k0,k1,lat) ! ! ps0 and embar0: do i=lon0,lon1 ps0(i,io2) = | b(i,1,1)*o2(lev0,i,lat)+b(i,1,2)*o1(lev0,i,lat)+fb(i,1) ps0(i,io1) = | b(i,2,1)*o2(lev0,i,lat)+b(i,2,2)*o1(lev0,i,lat)+fb(i,2) embar0(i) = 1./(ps0(i,io2)*rmassinv_o2+ps0(i,io1)*rmassinv_o1+ | (1.-ps0(i,io2)-ps0(i,io1))*rmassinv_n2) ! WKS4 = .5*(DMBAR/DZ)/MBAR wks4(i) = (embar(i,lev0)-embar0(i))/ | (dz*(embar0(i)+embar(i,lev0))) enddo ! i=lon0,lon1 ! write(6,"(/,'comp: lat=',i2,' ps0(:,io2)=',/,(6e12.4))") ! | lat,ps0(:,io2) ! write(6,"('ps0(:,io1)=',/,(6e12.4))") ps0(:,io1) ! write(6,"('embar0(:)=',/,(6e12.4))") embar0(:) ! write(6,"('wks4(:)=',/,(6e12.4))") wks4(:) ! ! ep, ak at level 1/2: km = 1 kp = 2 do i=lon0,lon1 ep(i,io2,kp) = 1.-(2./(embar0(i)+embar(i,lev0)))* | (rmass_o2+(embar(i,lev0)-embar0(i))/dz) ep(i,io1,kp) = 1.-(2./(embar0(i)+embar(i,lev0)))* | (rmass_o1+(embar(i,lev0)-embar0(i))/dz) zz(i,1,:) = 0. enddo ! i=lon0,lon1 ! write(6,"('ep(:,1,kp)=',/,(6e12.4))") ep(:,1,kp) ! write(6,"('ep(:,2,kp)=',/,(6e12.4))") ep(:,2,kp) do m=1,2 do i=lon0,lon1 ak(i,io2,m,kp) = | -delta(io2,m)*(phi(io1,3)+(phi(io1,io2)-phi(io1,3))* | .5*(ps0(i,io2)+o2(lev0,i,lat)))-(1.-delta(io2,m))* | (phi(io2,m)-phi(io2,3))*.5*(ps0(i,io2)+o2(lev0,i,lat)) ak(i,io1,m,kp) = | -delta(io1,m)*(phi(io2,3)+(phi(io2,io1)-phi(io2,3))* | .5*(ps0(i,io1)+o1(lev0,i,lat)))-(1.-delta(io1,m))* | (phi(io1,m)-phi(io1,3))*.5*(ps0(i,io1)+o1(lev0,i,lat)) enddo ! i=lon0,lon1 ! write(6,"('lat=',i2,' m=',i2,' ak(io2)=',/,(6e12.4))") ! | lat,m,ak(:,io2,m,kp) ! write(6,"('lat=',i2,' m=',i2,' ak(io1)=',/,(6e12.4))") ! | lat,m,ak(:,io1,m,kp) enddo ! m=1,2 ! ! WKS1=MBAR/M3*(T00/(T0+T))*0.25/(TAU*DET) AT LEVEL 1/2 ! tn lower boundary is stored in top slot tn(lev1..). do i=lon0,lon1 ! wks1(i) = 0.5*(embar0(i)+embar(i,lev0))*rmassinv_n2* ! | (t00/tn(lev1,i,lat))**0.25/(tau*(ak(i,1,1,kp)*ak(i,2,2,kp)- ! | ak(i,1,2,kp)*ak(i,2,1,kp))) ! Lower boundary of tn is now in tlbc: wks1(i) = 0.5*(embar0(i)+embar(i,lev0))*rmassinv_n2* | (t00/tlbc(i,lat))**0.25/(tau*(ak(i,1,1,kp)*ak(i,2,2,kp)- | ak(i,1,2,kp)*ak(i,2,1,kp))) enddo ! i=lon0,lon1 ! write(6,"('comp: lat=',i3,' wks1=',/,(6e12.4))") lat,wks1 ! ! Complete claculation of ak(1/2) do m=1,2 do i=lon0,lon1 ak(i,io2,m,kp) = ak(i,io2,m,kp)*wks1(i) ak(i,io1,m,kp) = ak(i,io1,m,kp)*wks1(i) gama(i,lev0,:,m) = 0. enddo ! i=lon0,lon1 ! write(6,"('lat=',i2,' m=',i2,' ak(io2)=',/,(6e12.4))") ! | lat,m,ak(:,io2,m,kp) ! write(6,"('lat=',i2,' m=',i2,' ak(io1)=',/,(6e12.4))") ! | lat,m,ak(:,io1,m,kp) enddo ! m=1,2 ! lonbeg = lon0 if (lon0==1) lonbeg = 3 lonend = lon1 if (lon1==nlonp4) lonend = nlonp4-2 ! ! Height (pressure) loop: ! For now (4/02), put k-loop on outside even tho embar and input ! fields are (k,i), for convenience in verification with tgcm15. ! km = 1 ! alternates 2,1,2,1,... during k-loop kp = 2 ! alternates 1,2,1,2,... during k-loop levloop: do k=lev0,lev1-1 ! DO 6 ktmp = km km = kp kp = ktmp do i=lon0,lon1 ep(i,io2,kp) = 1.-(2./(embar(i,k)+embar(i,k+1)))*(rmass_o2+ | (embar(i,k+1)-embar(i,k))/dz) ep(i,io1,kp) = 1.-(2./(embar(i,k)+embar(i,k+1)))*(rmass_o1+ | (embar(i,k+1)-embar(i,k))/dz) enddo ! i=lon0,lon1 ! write(6,"('lat=',i3,' k=',i3,' ep(io2)=',/,(6e12.4))") lat,k, ! | ep(:,io2,kp) ! write(6,"('ep(io1)=',/,(6e12.4))") ep(:,io1,kp) do m=1,2 do i=lon0,lon1 ! ! AK(K+1/2) ak(i,io2,m,kp) = | -delta(io2,m)*(phi(io1,3)+(phi(io1,io2)-phi(io1,3))* | .5*(o2(k,i,lat)+o2(k+1,i,lat)))- | (1.-delta(io2,m))*(phi(io2,m)-phi(io2,3))* | .5*(o2(k,i,lat)+o2(k+1,i,lat)) ak(i,io1,m,kp) = | -delta(io1,m)*(phi(io2,3)+(phi(io2,io1)-phi(io2,3))* | .5*(o1(k,i,lat)+o1(k+1,i,lat)))- | (1.-delta(io1,m))*(phi(io1,m)-phi(io1,3))* | .5*(o1(k,i,lat)+o1(k+1,i,lat)) enddo ! i=lon0,lon1 enddo ! m=1,2 ! ! WKS1=MBAR/M3*(T00/(T0+T))**0.25/(TAU*DET(ALFA)) do i=lon0,lon1 wks1(i) = 0.5*(embar(i,k)+embar(i,k+1))*rmassinv_n2* | (t00/(.5*(tn(k,i,lat)+tn(k+1,i,lat))))**0.25/ | (tau*(ak(i,1,1,kp)*ak(i,2,2,kp)-ak(i,1,2,kp)* | ak(i,2,1,kp))) ! ! EDDY DIFFUSION TERMS IN WKS3 AND WKS4 wks3(i) = wks4(i) wks4(i) = (embar(i,k+1)-embar(i,k))/ | (dz*(embar(i,k)+embar(i,k+1))) enddo ! i=lon0,lon1 ! ! FINISH CALCULATING AK(K+1/2) AND GENERATE PK, QK, RK do m=1,2 do isp=io2,io1 do i=lon0,lon1 ak(i,isp,m,kp) = ak(i,isp,m,kp)*wks1(i) pk(i,isp,m) = (ak(i,isp,m,km)*(1./dz+ep(i,m,km)/2.)- | expz(k)*(expzmid_inv*difk(k,iday)*dfactor(i)*(1./dz- | wks3(i))+0.25*(w(k,i,lat)+w(k+1,i,lat)))* | delta(isp,m))/dz rk(i,isp,m) = (ak(i,isp,m,kp)*(1./dz-ep(i,m,kp)/2.)- | expz(k)*(expzmid*difk(k+1,iday)*dfactor(i)*(1./dz+ | wks4(i))-0.25*(w(k,i,lat)+w(k+1,i,lat)))* | delta(isp,m))/dz qk(i,isp,m) = -(ak(i,isp,m,km)*(1./dz-ep(i,m,km)/2.)+ | ak(i,isp,m,kp)*(1./dz+ep(i,m,kp)/2.))/dz+expz(k)* | (((expzmid*difk(k+1,iday)*(1./dz-wks4(i))+expzmid_inv* | difk(k,iday)*(1./dz+wks3(i)))*dfactor(i)/dz+dtx2inv)* | delta(isp,m)-fs(i,k,isp,m,lat)) enddo ! i=lon0,lon1 ! write(6,"(/,'comp: lat=',i3,' m=',i2,' isp=',i2,' k=',i3)") ! | lat,m,isp,k ! write(6,"('ak=',/,(6e12.4))") ak(:,isp,m,kp) ! write(6,"('pk=',/,(6e12.4))") pk(:,isp,m) ! write(6,"('rk=',/,(6e12.4))") rk(:,isp,m) ! write(6,"('qk=',/,(6e12.4))") qk(:,isp,m) enddo ! isp=io2,io1 enddo ! m=1,2 ! ! Use advection saved from advecl calls at beginning of routine: ! do i=lon0,lon1 fk(i,io2) = o2_advec(k,i,lat) fk(i,io1) = o1_advec(k,i,lat) enddo ! i=lonbeg,lonend ! write(6,"(/,'comp: lat=',i3,' k=',i2)") lat,k ! write(6,"('fk(:,io2)=',/,(6e12.4))") fk(lonbeg:lonend,io2) ! write(6,"('fk(:,io1)=',/,(6e12.4))") fk(lonbeg:lonend,io1) ! ! Add explicit source terms to fk: do i=lon0,lon1 fk(i,io2) = fk(i,io2)-fs(i,k,io2,0,lat) fk(i,io1) = fk(i,io1)-fs(i,k,io1,0,lat) enddo ! i=lon0,lon1 ! write(6,"(/,'comp: lat=',i3,' k=',i2)") lat,k ! write(6,"('fk(:,io2)=',/,(6e12.4))") fk(lonbeg:lonend,io2) ! write(6,"('fk(:,io1)=',/,(6e12.4))") fk(lonbeg:lonend,io1) ! ! Complete calculation of rhs in fk: do i=lonbeg,lonend fk(i,io2) = expz(k)*(o2nm_smooth(k,i,lat)*dtx2inv-fk(i,io2)+ | hdo2(k,i,lat)) fk(i,io1) = expz(k)*(o1nm_smooth(k,i,lat)*dtx2inv-fk(i,io1)+ | hdo1(k,i,lat)) enddo ! i=lonbeg,lonend ! write(6,"(/,'comp: lat=',i3,' k=',i2)") lat,k ! write(6,"('fk(:,io2)=',/,(6e12.4))") fk(lonbeg:lonend,io2) ! write(6,"('fk(:,io1)=',/,(6e12.4))") fk(lonbeg:lonend,io1) ! ! fk is ok up to this point. ! In earlier version, periodic points for fk were taken here. ! For now, ignore periodic points. ! ! Lower boundary: if (k==lev0) then do m=1,2 ! DO 16 do kk=1,2 do i=lon0,lon1 qk(i,io2,m) = qk(i,io2,m)+pk(i,io2,kk)*b(i,kk,m) qk(i,io1,m) = qk(i,io1,m)+pk(i,io1,kk)*b(i,kk,m) enddo ! i=lon0,lon1 enddo ! kk=1,2 enddo ! m=1,2 do m=1,2 do i=lon0,lon1 fk(i,io2) = fk(i,io2)-pk(i,io2,m)*fb(i,m) fk(i,io1) = fk(i,io1)-pk(i,io1,m)*fb(i,m) pk(i,:,m) = 0. enddo ! i=lon0,lon1 enddo ! m=1,2 ! do m=1,2 ! write(6,"('comp lbc: m=',i2,' lat=',i2)") m,lat ! write(6,"('qk(io2)=',/,(6e12.4))") qk(:,io2,m) ! write(6,"('qk(io1)=',/,(6e12.4))") qk(:,io1,m) ! enddo ! m=1,2 ! write(6,"('fk(io2)=',/,(6e12.4))") fk(:,io2) ! write(6,"('fk(io1)=',/,(6e12.4))") fk(:,io1) ! ! Upper boundary: elseif (k==lev1-1) then do m=1,2 do i=lon0,lon1 qk(i,io2,m) = qk(i,io2,m)+(1.+.5*ep(i,m,kp)*dz)/ | (1.-.5*ep(i,m,kp)*dz)*rk(i,io2,m) qk(i,io1,m) = qk(i,io1,m)+(1.+.5*ep(i,m,kp)*dz)/ | (1.-.5*ep(i,m,kp)*dz)*rk(i,io1,m) rk(i,:,m) = 0. enddo ! i=lon0,lon1 enddo ! m=1,2 ! do m=1,2 ! write(6,"('comp ubc: m=',i2,' lat=',i2)") m,lat ! write(6,"('qk(io2)=',/,(6e12.4))") qk(:,io2,m) ! write(6,"('qk(io1)=',/,(6e12.4))") qk(:,io1,m) ! enddo ! m=1,2 endif ! lbc or ubc ! ! QK=ALFAK=QK-PK*GAMA(K-1) do m=1,2 ! DO 18 do kk=1,2 do i=lon0,lon1 ! write(6,"('comp: i=',i2,' kk=',i2,' m=',i2,' k=',i2, ! | ' lat=',i2)") i,kk,m,k,lat ! write(6,"('qk=',e12.4,' pk=',e12.4,' gama=',e12.4))") ! | qk(i,io2,m),pk(i,io2,kk),gama(i,k,kk,m) qk(i,io2,m) = qk(i,io2,m)-pk(i,io2,kk)*gama(i,k,kk,m) qk(i,io1,m) = qk(i,io1,m)-pk(i,io1,kk)*gama(i,k,kk,m) enddo ! i=lon0,lon1 enddo ! kk=1,2 enddo ! m=1,2 ! Testing will not work here until gama is incremented below.. ! do m=1,2 ! write(6,"('comp: m=',i2,' k=',i2,' lat=',i2)") m,k,lat ! write(6,"('qk(io2)=',/,(6e12.4))") qk(:,io2,m) ! write(6,"('qk(io1)=',/,(6e12.4))") qk(:,io1,m) ! enddo ! m=1,2 ! ! WKS1=DET(ALFA) do i=lon0,lon1 wks1(i) = qk(i,1,1)*qk(i,2,2)-qk(i,1,2)*qk(i,2,1) enddo ! i=lon0,lon1 ! ! WKM1=ALFAI do m=1,2 do i=lon0,lon1 wkm1(i,io2,m) = (delta(io2,m)*qk(i,io1,io1)- | (1.-delta(io2,m))*qk(i,io2,m))/wks1(i) wkm1(i,io1,m) = (delta(io1,m)*qk(i,io2,io2)- | (1.-delta(io1,m))*qk(i,io1,m))/wks1(i) enddo ! i=lon0,lon1 enddo ! m=1,2 ! ! WKV1=FK-PK*Z(K) do i=lon0,lon1 wkv1(i,io2) = fk(i,io2) wkv1(i,io1) = fk(i,io1) enddo ! i=lon0,lon1 ! ! GAMA(K+1)=ALFAI*RK do m=1,2 do i=lon0,lon1 gama(i,k+1,io2,m) = 0. gama(i,k+1,io1,m) = 0. wkv1(i,io2) = wkv1(i,io2)-pk(i,io2,m)*zz(i,k,m) wkv1(i,io1) = wkv1(i,io1)-pk(i,io1,m)*zz(i,k,m) enddo ! i=lon0,lon1 do kk=1,2 do i=lon0,lon1 gama(i,k+1,io2,m) = gama(i,k+1,io2,m)+wkm1(i,io2,kk)* | rk(i,kk,m) gama(i,k+1,io1,m) = gama(i,k+1,io1,m)+wkm1(i,io1,kk)* | rk(i,kk,m) enddo ! i=lon0,lon1 enddo ! kk=1,2 enddo ! m=1,2 ! ! Z(K+1)=WKM1*WKV1 do i=lon0,lon1 zz(i,k+1,:) = 0. enddo ! i=lon0,lon1 do m=1,2 do i=lon0,lon1 zz(i,k+1,io2) = zz(i,k+1,io2)+wkm1(i,io2,m)*wkv1(i,m) zz(i,k+1,io1) = zz(i,k+1,io1)+wkm1(i,io1,m)*wkv1(i,m) enddo ! i=lon0,lon1 enddo ! m=1,2 ! do m=1,2 ! write(6,"('comp: m=',i2,' k=',i2,' lat=',i2)") m,k,lat ! write(6,"('gama(k+1,io2)=',/,(6e12.4))") gama(:,k+1,io2,m) ! write(6,"('gama(k+1,io1)=',/,(6e12.4))") gama(:,k+1,io1,m) ! enddo ! write(6,"('zz(k+1,io2)=',/,(6e12.4))") zz(:,k+1,io2) ! write(6,"('zz(k+1,io1)=',/,(6e12.4))") zz(:,k+1,io1) ! ! End main pressure loop: enddo levloop ! k=lev0,lev1-1 ! ! Save diagnostics: ! real,dimension(lev0:lev1,lon0:lon1,2) :: zz_ki ! real,dimension(lev0:lev1,lon0:lon1,2,2) :: gama_ki ! do k=lev0,lev1 zz_ki(k,:,io2) = zz(:,k,io2) zz_ki(k,:,io1) = zz(:,k,io1) do m=1,2 gama_ki(k,:,io2,m) = gama(:,k,io2,m) gama_ki(k,:,io1,m) = gama(:,k,io1,m) enddo enddo ! call addfld('ZZ_O2',' ',' ',zz_ki(k0:k1-1,:,io2), ! | 'lev',k0,k1-1,'lon',i0,i1,lat) ! call addfld('ZZ_O1',' ',' ',zz_ki(k0:k1-1,:,io1), ! | 'lev',k0,k1-1,'lon',i0,i1,lat) ! call addfld('GAMAO2M1',' ',' ',gama_ki(k0:k1-1,:,io2,1), ! | 'lev',k0,k1-1,'lon',i0,i1,lat) ! call addfld('GAMAO2M2',' ',' ',gama_ki(k0:k1-1,:,io2,2), ! | 'lev',k0,k1-1,'lon',i0,i1,lat) ! call addfld('GAMAO1M1',' ',' ',gama_ki(k0:k1-1,:,io1,1), ! | 'lev',k0,k1-1,'lon',i0,i1,lat) ! call addfld('GAMAO1M2',' ',' ',gama_ki(k0:k1-1,:,io1,2), ! | 'lev',k0,k1-1,'lon',i0,i1,lat) ! ! Set upper boundary to zero: do i=lon0,lon1 o2_upd(lev1,i,lat) = 0. o1_upd(lev1,i,lat) = 0. upd(lev1,i,:) = 0. enddo ! i=lon0,lon1 ! ! Downward sweep: do k=lev1-1,lev0,-1 do isp=io2,io1 do i=lon0,lon1 upd(k,i,isp) = zz(i,k+1,isp) enddo ! i=lon0,lon1 do m=1,2 do i=lon0,lon1 upd(k,i,isp) = upd(k,i,isp)-gama(i,k+1,isp,m)* | upd(k+1,i,m) enddo enddo ! m=1,2 enddo ! isp=io2,io1 enddo ! k=lev1-1,lev0,-1 ! ! Transfer to output arrays: do k=lev0,lev1 o2_upd(k,lon0:lon1,lat) = upd(k,:,io2) o1_upd(k,lon0:lon1,lat) = upd(k,:,io1) enddo ! k=lev0,lev1 ! ! Upper boundary: ! kp is carried forward from the last iteration of levloop above. do i=lon0,lon1 o2_upd(lev1,i,lat) = | (1.+.5*ep(i,io2,kp)*dz)/ | (1.-.5*ep(i,io2,kp)*dz)*o2_upd(lev1-1,i,lat) o1_upd(lev1,i,lat) = | (1.+.5*ep(i,io1,kp)*dz)/ | (1.-.5*ep(i,io1,kp)*dz)*o1_upd(lev1-1,i,lat) enddo ! i=lon0,lon1 ! ! call addfld('O2_SOLV',' ',' ',o2_upd(:,lon0:lon1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! call addfld('O1_SOLV',' ',' ',o1_upd(:,lon0:lon1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! ! End latitude scan: enddo ! lat=lat0,lat1 ! ! Filter the new composition species: ! ! Fourier smoothing of O2 and O: call filter_o2o(o2_upd,lev0,lev1,lon0,lon1,lat0,lat1) call filter_o2o(o1_upd,lev0,lev1,lon0,lon1,lat0,lat1) ! ! Resume latitude scan: do lat=lat0,lat1 ! call addfld('O2_FILT',' ',' ',o2_upd(:,lon0:lon1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! call addfld('O1_FILT',' ',' ',o1_upd(:,lon0:lon1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! ! Time smoothing: do i=lon0,lon1 do k=lev0,lev1 o2nm_upd(k,i,lat) = dtsmooth_div2*(o2_nm(k,i,lat)+ | o2_upd(k,i,lat)) + dtsmooth*o2(k,i,lat) o1nm_upd(k,i,lat) = dtsmooth_div2*(o1_nm(k,i,lat)+ | o1_upd(k,i,lat)) + dtsmooth*o1(k,i,lat) enddo ! k=lev0,lev1-1 enddo ! i=lon0,lon1 ! call addfld('O2NM_OUT',' ',' ',o2nm_upd(:,lon0:lon1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! call addfld('O1NM_OUT',' ',' ',o1nm_upd(:,lon0:lon1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) enddo ! lat=lat0,lat1 #ifdef MPI ! ! Periodic points: ! call mp_periodic_f3d(o2_upd(:,lon0:lon1,lat0-1:lat1+1), ! | lev0,lev1,lon0,lon1,lat0-1,lat1+1) ! call mp_periodic_f3d(o1_upd(:,lon0:lon1,lat0-1:lat1+1), ! | lev0,lev1,lon0,lon1,lat0-1,lat1+1) #endif ! ! Insure non-negative O2,O: do lat=lat0,lat1 do i=lon0,lon1 do k=lev0,lev1 if (o2_upd(k,i,lat) < small) o2_upd(k,i,lat) = small if (o1_upd(k,i,lat) < small) o1_upd(k,i,lat) = small if (o2nm_upd(k,i,lat) < small) o2nm_upd(k,i,lat) = small if (o1nm_upd(k,i,lat) < small) o1nm_upd(k,i,lat) = small if (1.-small-o2_upd(k,i,lat)-o1_upd(k,i,lat) < 0.) then o2_upd(k,i,lat) = o2_upd(k,i,lat)*((1.-small)/ | (o2_upd(k,i,lat)+o1_upd(k,i,lat))) o1_upd(k,i,lat) = o1_upd(k,i,lat)*((1.-small)/ | (o2_upd(k,i,lat)+o1_upd(k,i,lat))) endif if (1.-small-o2nm_upd(k,i,lat)-o1nm_upd(k,i,lat) < 0.) then o2nm_upd(k,i,lat) = o2nm_upd(k,i,lat)*((1.-small)/ | (o2nm_upd(k,i,lat)+o1nm_upd(k,i,lat))) o1nm_upd(k,i,lat) = o1nm_upd(k,i,lat)*((1.-small)/ | (o2nm_upd(k,i,lat)+o1nm_upd(k,i,lat))) endif enddo ! k=lev0,lev1 enddo ! i=lon0,lon1 call mkdiag_O_N2('O_N2',o1_upd(:,lon0:lon1,lat), | o2_upd(:,lon0:lon1,lat),lev0,lev1,lon0,lon1,lat) ! call addfld('O2_OUT',' ',' ',o2_upd(:,i0:i1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! call addfld('O1_OUT',' ',' ',o1_upd(:,i0:i1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! call addfld('O2_NMOUT',' ',' ',o2nm_upd(:,i0:i1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! call addfld('O1_NMOUT',' ',' ',o1nm_upd(:,i0:i1,lat), ! | 'lev',k0,k1,'lon',i0,i1,lat) ! enddo ! lat=lat0,lat1 end subroutine comp !----------------------------------------------------------------------- subroutine advecl(o2,o1,un,vn,o2_advec,o1_advec, | lev0,lev1,lon0,lon1,lat0,lat1,lat) ! ! Horizontal advection for O2,O. ! In previous versions, this was sub advecl (inline.F), called from ! k-loop in comp.F. Here it is called from beginning of comp.F at ! all latitudes and saved for later use (fk) inside comp.F k-loop. ! O2,o1,un,vn already have i-2,i-1,i+1,i+2, and j-1,j-2,j+1,j+2 for ! finite differencing. ! use cons_module,only: dlamda_2div3 ,dlamda_1div12, dphi_2div3, | dphi_1div12,re_inv,racs implicit none ! ! Args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1,lat real,dimension(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2),intent(in):: | o2,o1,un,vn real,dimension(lev0:lev1,lon0:lon1,lat0:lat1),intent(out) :: | o2_advec, o1_advec ! ! Local: integer :: k,i ! do k=lev0,lev1-1 do i=lon0,lon1 o2_advec(k,i,lat) = .5*racs(lat)* | (dlamda_2div3*(o2(k,i+1,lat)-o2(k,i-1,lat))* | (un(k,i+1,lat)+un(k,i-1,lat))- | dlamda_1div12*(o2(k,i+2,lat)-o2(k,i-2,lat))* | (un(k,i+2,lat)+un(k,i-2,lat)))+ | .5*re_inv* | (dphi_2div3*(o2(k,i,lat+1)-o2(k,i,lat-1))* | (vn(k,i,lat+1)+vn(k,i,lat-1))- | dphi_1div12*(o2(k,i,lat+2)-o2(k,i,lat-2))* | (vn(k,i,lat+2)+vn(k,i,lat-2))) o1_advec(k,i,lat) = .5*racs(lat)* | (dlamda_2div3*(o1(k,i+1,lat)-o1(k,i-1,lat))* | (un(k,i+1,lat)+un(k,i-1,lat))- | dlamda_1div12*(o1(k,i+2,lat)-o1(k,i-2,lat))* | (un(k,i+2,lat)+un(k,i-2,lat)))+ | .5*re_inv* | (dphi_2div3*(o1(k,i,lat+1)-o1(k,i,lat-1))* | (vn(k,i,lat+1)+vn(k,i,lat-1))- | dphi_1div12*(o1(k,i,lat+2)-o1(k,i,lat-2))* | (vn(k,i,lat+2)+vn(k,i,lat-2))) enddo ! i=lon0,lon1 enddo ! k=lev0,lev1-1 end subroutine advecl !----------------------------------------------------------------------- subroutine filter_o2o(fout,lev0,lev1,lon0,lon1,lat0,lat1) ! ! Filter updated W omega: ! use params_module,only: nlat,nlonp4,nlon use filter_module,only: filter use cons_module,only: kut ! kut(nlat) #ifdef MPI use mpi_module,only: mp_gatherlons_f3d,mp_scatterlons_f3d,mytidi implicit none #else implicit none integer :: mytidi=0 #endif ! ! Args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1 real,intent(inout) :: fout(lev0:lev1,lon0-2:lon1+2,lat0-2:lat1+2) ! ! VT vampir tracing: ! #ifdef VT #include #endif ! ! Local: integer :: i,j,nlevs,nlons,nlats real :: fik(nlonp4,lev0:lev1),fkij(lev0:lev1,nlonp4,lat0:lat1) real :: fmin,fmax ! #ifdef VT ! code = 131 ; state = 'filter_o2o' ; activity='Filtering' call vtbegin(131,ier) #endif ! nlevs = lev1-lev0+1 nlons = lon1-lon0+1 nlats = lat1-lat0+1 ! ! Define lons in w_ki from current task: fkij = 0. do j=lat0,lat1 do i=lon0,lon1 fkij(:,i,j) = fout(:,i,j) enddo enddo ! j=lat0,lat1 ! #ifdef MPI ! ! Gather longitudes into tasks in first longitude column of task table ! (leftmost of each j-row) for global fft. (i.e., tasks with mytidi==0 ! gather lons from other tasks in that row). This includes all latitudes. ! call mp_gatherlons_f3d(fkij,lev0,lev1,lon0,lon1,lat0,lat1,1) #endif ! ! Only leftmost tasks at each j-row of tasks does the global filtering: if (mytidi==0) then ! ! Define 2d array with all longitudes for filter at each latitude: latscan: do j=lat0,lat1 do i=1,nlonp4 fik(i,:) = fkij(:,i,j) enddo ! i=1,nlonp4 ! ! Remove wave numbers > kut(lat): call filter(fik,lev0,lev1,kut(j),j) ! tiegcm ! call filter2(fik,lev0,lev1,j) ! timegcm ! ! Return filtered array to fkij: do i=1,nlonp4 fkij(:,i,j) = fik(i,:) enddo ! i=1,nlonp4 enddo latscan ! j=lat0,lat1 endif ! mytidi==0 #ifdef MPI ! ! Now leftmost task at each j-row must redistribute filtered data ! back to other tasks in the j-row (mytidi>0,mytidj) (includes latitude): ! call mp_scatterlons_f3d(fkij,lev0,lev1,lon0,lon1,lat0,lat1,1) #endif ! ! Return filtered array to fout at current task longitudes and latitudes: do j=lat0,lat1 do i=lon0,lon1 fout(:,i,j) = fkij(:,i,j) enddo enddo ! #ifdef VT ! code = 131 ; state = 'filter_o2o' ; activity='Filtering' call vtend(131,ier) #endif end subroutine filter_o2o