SUBROUTINE SW_venus_dc( PRMU0, PFRAC, S PPB, pt, S PHEAT, S PTOPSW,PSOLSW,ZFSNET) use dimphy IMPLICIT none #include "dimensions.h" #include "YOMCST.h" C C ------------------------------------------------------------------ C C PURPOSE. C -------- C c this routine loads and interpolates the shortwave radiation c fluxes taken from Dave Crisp calculations for Venus. c Ref: Crisp 1986. C C AUTHOR. C ------- C Sebastien Lebonnois C C MODIFICATIONS. C -------------- C ORIGINAL : 27/07/2005 C ------------------------------------------------------------------ C C* ARGUMENTS: C c inputs REAL PRMU0 ! COSINE OF ZENITHAL ANGLE REAL PFRAC ! fraction de la journee REAL PPB(KFLEV+1) ! inter-couches PRESSURE (bar) REAL pt(kflev) ! mid-layer temperature C c output REAL PHEAT(KFLEV) ! SHORTWAVE HEATING (K/VENUSDAY) within each layer REAL PTOPSW ! SHORTWAVE FLUX AT T.O.A. (net) REAL PSOLSW ! SHORTWAVE FLUX AT SURFACE (net) REAL ZFSNET(KFLEV+1) ! net solar flux at ppb levels C C* LOCAL VARIABLES: C integer nldc,nszadc real dureejour parameter (nldc=49) ! fichiers Crisp parameter (nszadc=8) ! fichiers Crisp parameter (dureejour=10.087e6) integer i,j,nsza,nsza0,nl0 real solarrate ! solar heating rate (K/earthday) real zsnet(nldc+1,nszadc) ! net solar flux (W/m**2) (+ vers bas) real zsdn,zsup ! downward/upward solar flux (W/m**2) real solza(nszadc) ! solar zenith angles in table real presdc(nldc+1) ! pressure levels in table (bar) real tempdc(nldc+1) ! temperature in table (K) real altdc(nldc+1) ! altitude in table (km) real coolrate ! IR heating rate (K/earthday) ? real totalrate ! total rate (K/earthday) real zldn ! downward IR flux (W/m**2) ? real zlup ! upward IR flux (W/m**2) ? real zsolnet(nldc) ! for testing mean net solar flux in DC character*22 nullchar real sza0,factsza,factflux logical firstcall data firstcall/.true./ REAL,save,allocatable :: zsolVE(:) ! net solar flux at ppb levels, fichiers VE save solza,zsnet,presdc,tempdc,altdc save firstcall c ------------------------ c Loading the file c ------------------------ if (firstcall) then allocate(zsolVE(klevp1)) open(11,file='dataDCrisp.dat') read(11,*) nullchar do nsza=1,nszadc read(11,*) nullchar read(11,*) nullchar read(11,*) nullchar read(11,'(22x,F11.5)') solza(nsza) read(11,*) nullchar read(11,*) nullchar read(11,*) nullchar read(11,'(3(2x,F10.4),36x,4(2x,F11.5))') . presdc(nldc+1),tempdc(nldc+1), altdc(nldc+1), . zsdn,zsup,zldn,zlup zsnet(nldc+1,nsza)=zsdn-zsup do i=1,nldc j = nldc+1-i ! changing: vectors from surface to top read(11,'(6(2x,F10.4),4(2x,F11.5))') . presdc(j),tempdc(j),altdc(j), . solarrate,coolrate,totalrate, . zsdn,zsup,zldn,zlup zsnet(j,nsza)=zsdn-zsup enddo enddo close(11) c ----------- TEST ------------ c Fichiers de Vincent c ----------------------------- c open(12,file='flux_vis_dcGCM.txt') c read(12,*) nullchar c c do j=1,kflev+1 c read(12,*) zlup,zldn,zsolVE(j) c enddo c c close(12) c ----------------------------- c -------- FIN TEST ---------- firstcall=.false. endif c ----------- TEST ------------ c Moyenne planetaire c ----------------------------- c do j=1,nldc c --- c zsolnet(j) = zsnet(j,1)*0.5* c . sin(solza(1)*RPI/180.)*solza(2)*RPI/180./2. c do nsza=2,nszadc-1 c zsolnet(j) = zsolnet(j)+zsnet(j,nsza)*0.5* c . sin(solza(nsza)*RPI/180.)* c . (solza(nsza+1)-solza(nsza-1))*RPI/180./2. c enddo c zsolnet(j) = zsolnet(j)+zsdn(j,nszadc)*0.5* c . sin(solza(nszadc)*RPI/180.)* c . (90.-solza(nszadc-1))*RPI/180./2. c --- c zsolnet(j) = 0.0 c do nsza=1,nszadc-1 c zsolnet(j) = zsolnet(j)+(zsnet(j,nsza ) c . +zsnet(j,nsza+1))*0.5* c . (cos(solza(nsza )*RPI/180.)- c . cos(solza(nsza+1)*RPI/180.) ) c enddo c zsolnet(j) = zsolnet(j)+zsnet(j,nszadc)*0.25* c . cos(solza(nszadc)*RPI/180.) c --- c print*,j,altdc(j),zsolnet(j) c enddo c stop c ----------------------------- c -------- FIN TEST ---------- c -------------------------------------- c Interpolation in the GCM vertical grid c -------------------------------------- c Zenith angle c ------------ sza0 = acos(PRMU0)/3.1416*180. c print*,'Angle Zenithal =',sza0,' PFRAC=',PFRAC do nsza=1,nszadc if (solza(nsza).le.sza0) then nsza0 = nsza+1 endif enddo if (nsza0.ne.nszadc+1) then factsza = (sza0-solza(nsza0-1))/(solza(nsza0)-solza(nsza0-1)) else factsza = min((sza0-solza(nszadc))/(90.-solza(nszadc)), 1.) endif c Pressure levels c --------------- do j=1,kflev+1 nl0 = 2 do i=1,nldc if (presdc(i).ge.PPB(j)) then nl0 = i+1 endif enddo factflux = (log10(max(PPB(j),presdc(nldc+1))) . -log10(presdc(nl0-1))) . /(log10(presdc(nl0))-log10(presdc(nl0-1))) c factflux = (max(PPB(j),presdc(nldc+1))-presdc(nl0-1)) c . /(presdc(nl0)-presdc(nl0-1)) if (nsza0.ne.nszadc+1) then ZFSNET(j) = factflux * factsza *zsnet(nl0,nsza0) . + factflux *(1.-factsza)*zsnet(nl0,nsza0-1) . + (1.-factflux)* factsza *zsnet(nl0-1,nsza0) . + (1.-factflux)*(1.-factsza)*zsnet(nl0-1,nsza0-1) else ZFSNET(j) = factflux *(1.-factsza)*zsnet(nl0,nsza0-1) . + (1.-factflux)*(1.-factsza)*zsnet(nl0-1,nsza0-1) endif ZFSNET(j) = ZFSNET(j)*PFRAC enddo c ----------- TEST ------------ c Fichiers de Vincent c ----------------------------- c do j=1,kflev+1 c ZFSNET(j)=zsolVE(j) c enddo c ----------------------------- c -------- FIN TEST ---------- PTOPSW = ZFSNET(kflev+1) PSOLSW = ZFSNET(1) c Heating rates c ------------- c On utilise le gradient du flux pour calculer le taux de chauffage: c heat(K/s) = d(fluxnet) (W/m2) c *g (m/s2) c /(-dp) (epaisseur couche, en Pa=kg/m/s2) c /cp (J/kg/K) do j=1,kflev ! ADAPTATION GCM POUR CP(T) PHEAT(j) = (ZFSNET(j+1)-ZFSNET(j)) . *RG/cpdet(pt(j)) / ((PPB(j)-PPB(j+1))*1.e5) PHEAT(j) = PHEAT(j)*dureejour ! K/venus_day enddo return end