source: trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh.F @ 3567

Last change on this file since 3567 was 2503, checked in by emillour, 4 years ago

Venus GCM:
Small correction in the computation of SW fluxes for the
very specific case when running without a diurnal cycle.
VB+EM
File size: 10.2 KB
RevLine 
[1591]1      SUBROUTINE SW_venus_rh(PRMU0, PFRAC, latdeg,
[1723]2     S              PPA, PPB, pt,
[1591]3     S              PHEAT,
4     S              PTOPSW,PSOLSW,ZFSNET)
5     
[2503]6      use dimphy, only: klev
[1621]7      use cpdet_phy_mod, only: cpdet
[1591]8      IMPLICIT none
9
[2503]10      include "YOMCST.h"
11      include "clesphys.h"
[1591]12C
13C     ------------------------------------------------------------------
14C
15C     PURPOSE.
16C     --------
17C
18c      this routine loads and interpolates the shortwave radiation
19c     fluxes taken from Rainer Haus calculations for Venus.
20c     Ref: Haus et al. 2016
21C
22C     AUTHOR.
23C     -------
24C        Sebastien Lebonnois
25C
26C     MODIFICATIONS.
27C     --------------
28C        ORIGINAL : 5/2016
29C     ------------------------------------------------------------------
30C
31C* ARGUMENTS:
32C
33c inputs
34
[2503]35      REAL,INTENT(IN) ::   PRMU0  ! COSINE OF ZENITHAL ANGLE
36      REAL,INTENT(IN) ::   PFRAC  ! fraction de la journee
37      REAL,INTENT(IN) ::   latdeg ! |latitude| (in degrees)
38      REAL,INTENT(IN) ::   PPB(klev+1)  ! inter-couches PRESSURE (bar)
39      REAL,INTENT(IN) ::   PPA(klev)
40      REAL,INTENT(IN) ::   pt(klev)     ! mid-layer temperature
[1591]41C
42c output
43
[2503]44      REAL,INTENT(OUT) ::   PHEAT(klev)  ! SHORTWAVE HEATING (K/s) within each layer
[1723]45      REAL   PHEATPPA(klev)
[2503]46      REAL,INTENT(OUT) ::   PTOPSW       ! SHORTWAVE FLUX AT T.O.A. (net)
47      REAL,INTENT(OUT) ::   PSOLSW       ! SHORTWAVE FLUX AT SURFACE (net)
48      REAL,INTENT(OUT) ::   ZFSNET(klev+1) ! net solar flux at ppb levels
[1591]49
50C
51C* LOCAL VARIABLES:
52C
53      integer nlrh,nszarh,nlatrh
54      parameter (nlrh=118)  ! fichiers Rainer Haus
55      parameter (nszarh=7) ! fichiers Rainer Haus
56      parameter (nlatrh=19) ! fichiers Rainer Haus
57     
[1723]58      integer i,j,k,lat,nsza,nsza0(2),nl0,nlat0
[1591]59      real   zsnet(nlrh+1,nszarh+1,nlatrh+1)! net solar flux (W/m**2) (+ vers bas)
60      real   solza(nszarh,nlatrh)       ! solar zenith angles in table
61      real   presrh(nlrh+1)             ! pressure in table (bar)
[1726]62      real   logplayrh(nlrh)
[1591]63      real   altrh(nlrh+1)              ! altitude in table (km)
64      real   latrh(nlatrh)              ! latitude in table (degrees)
65      character*22 nullchar
66      real   sza0,factsza(2),factflux,factlat
67      real   zsnetmoy
68      logical firstcall
69      data    firstcall/.true./
70      save   solza,zsnet,altrh,latrh,presrh
71      save   firstcall
[1723]72      real   Tplay(nlrh)
[1726]73      real   Qrh1(nlrh)
74      real   Qrh2(nlrh)
75      real   Qrh3(nlrh)
76      real   Qrh4(nlrh)
[1591]77     
78c ------------------------
79c Loading the file
80c ------------------------
81      if (firstcall) then
82
83       zsnet=0.
84
85       open(11,file='SolarNetFlux_RH.dat')
86
87       do i=1,nlrh+1
88          read(11,'(E5.1,4x,F8.2)') altrh(i),presrh(i)
89       enddo
90
91       do lat=1,nlatrh
92         latrh(lat)=5.*(lat-1)
93         read(11,*) nullchar
94         read(11,*) nullchar
95         read(11,'(3x,7(5x,E8.5))') solza(:,lat)
96         read(11,*) nullchar
97
98         do i=1,nlrh+1
99          read(11,'(E6.1,7(2x,F11.5),7x,F11.5)')
100     .          altrh(i),zsnet(i,1:nszarh,lat),zsnetmoy
101         enddo
102         read(11,*) nullchar
103       enddo
104       latrh(nlatrh)=89.
105
106c Correction of factor 2 in the table...
107       zsnet=zsnet*2.
108
109       close(11)
110
111       firstcall=.false.
112      endif
113
114c --------------------------------------
115c Interpolation in the GCM vertical grid
116c --------------------------------------
117
118c Latitude
119c ---------
120     
121      do lat=1,nlatrh
122         if (latrh(lat).le.latdeg) then
123              nlat0 = lat+1
124         endif
125      enddo
[1723]126
[1591]127      if (nlat0.ne.nlatrh+1) then
128        factlat = (latdeg-latrh(nlat0-1))/(latrh(nlat0)-latrh(nlat0-1))
129      else
130        factlat = min((latdeg-latrh(nlatrh))/(90.-latrh(nlatrh)), 1.)
131      endif
132
133c Zenith angle
134c ------------
135     
136      sza0 = acos(PRMU0)/3.1416*180.
137      nsza0(:)=2
138
[2503]139      if (.not.cycle_diurne) then
140        ! without a diurnal cycle, no need for any elaborate weights of sza
141        factsza(1)=1
142        factsza(2)=0
143      else
144        ! standard case with diurnal cycle
145        do nsza=1,nszarh
146          if (solza(nsza,nlat0-1).le.sza0) then
[1591]147              nsza0(1) = nsza+1
[2503]148          endif
149        enddo
150        if (nsza0(1).ne.nszarh+1) then
151            factsza(1) = (sza0-solza(nsza0(1)-1,nlat0-1))/
152     .          (solza(nsza0(1),nlat0-1)-solza(nsza0(1)-1,nlat0-1))
153        else
154            factsza(1) = min((sza0-solza(nszarh,nlat0-1))/
155     .           (90.-solza(nszarh,nlat0-1)), 1.)
156        endif
157        if (nlat0.ne.nlatrh+1) then
158          do nsza=1,nszarh
159            if (solza(nsza,nlat0).le.sza0) then
[1591]160              nsza0(2) = nsza+1
[2503]161            endif
162          enddo
163          if (nsza0(2).eq.nszarh+1) then
164             factsza(2) = min((sza0-solza(nszarh,nlat0))/
165     .           (90.-solza(nszarh,nlat0)), 1.)
166          elseif ((nsza0(2).eq.2).and.(solza(1,nlat0).gt.sza0)) then
167            factsza(2) = 0.
168          else
169            factsza(2) = (sza0-solza(nsza0(2)-1,nlat0))/
170     .          (solza(nsza0(2),nlat0)-solza(nsza0(2)-1,nlat0))
171          endif
172        else
173          nsza0(2)   = nszarh+1
174          factsza(2) = 1.
175        endif ! of if (nlat0.ne.nlatrh+1)
176      endif ! of if (.not.cycle_diurne)
177
[1591]178c Pressure levels
179c ---------------
180      do j=1,klev+1
181        nl0 = nlrh
182        do i=nlrh+1,2,-1
183           if (presrh(i).ge.PPB(j)) then
184                nl0 = i-1
185           endif
186        enddo
187
188        factflux = (log10(max(PPB(j),presrh(1)))-log10(presrh(nl0+1)))
189     .            /(log10(presrh(nl0))-log10(presrh(nl0+1)))
190
191        ZFSNET(j) =  factlat*(
192     .      factflux   *  factsza(2)   *zsnet(nl0,nsza0(2),nlat0)
193     . +   factflux   *(1.-factsza(2))*zsnet(nl0,nsza0(2)-1,nlat0)
194     . + (1.-factflux)*  factsza(2)   *zsnet(nl0+1,nsza0(2),nlat0)
195     . + (1.-factflux)*(1.-factsza(2))*zsnet(nl0+1,nsza0(2)-1,nlat0) )
196     .            + (1.-factlat)*(
197     .      factflux   *  factsza(1)   *zsnet(nl0,nsza0(1),nlat0-1)
198     . +   factflux   *(1.-factsza(1))*zsnet(nl0,nsza0(1)-1,nlat0-1)
199     . + (1.-factflux)*  factsza(1)   *zsnet(nl0+1,nsza0(1),nlat0-1)
200     . + (1.-factflux)*(1.-factsza(1))*zsnet(nl0+1,nsza0(1)-1,nlat0-1) )
201
202        ZFSNET(j) = ZFSNET(j)*PFRAC
203       
204      enddo
205      PTOPSW = ZFSNET(klev+1)
206      PSOLSW = ZFSNET(1)
[1723]207
208#ifdef MESOSCALE
[1726]209! extrapolation play RH pressure
[1723]210      do j=1,nlrh
[1726]211         logplayrh(j)=(log(presrh(j+1))+log(presrh(j)))/2.
[1723]212      enddo
[1726]213! Extrapolation of temperature over RH play pressure
[1723]214      do i=nlrh,2,-1
215        nl0 = 2
216        do j=1,klev-1
[1726]217           if (exp(logplayrh(i)).le.PPA(j)) then
[1723]218                nl0 = j+1
219           endif
220        enddo
[1726]221        factflux = (log10(max(exp(logplayrh(i)),PPA(klev)))
[1723]222     .                         -log10(PPA(nl0-1)))
223     .       /(log10(PPA(nl0))-log10(PPA(nl0-1)))
224        Tplay(i)=factflux*pt(nl0)
225     .             + (1.-factflux)*pt(nl0-1)
226
227      ENDDO
[1726]228! RH PHEAT over RH play pressure
[1723]229      DO k=1,nlrh
230c
[1726]231       Qrh1(k)=((RG/cpdet(Tplay(k)))
[1723]232     .     *((zsnet(k+1,nsza0(1),nlat0-1)-zsnet(k,nsza0(1),nlat0-1))
233     .         *PFRAC))
234     .      /((presrh(k)-presrh(k+1))*1.e5)
[1726]235       Qrh2(k)=((RG/cpdet(Tplay(k)))
[1723]236     . *((zsnet(k+1,nsza0(1)-1,nlat0-1)-zsnet(k,nsza0(1)-1,nlat0-1))
237     .         *PFRAC))
238     .      /((presrh(k)-presrh(k+1))*1.e5)
[1726]239       Qrh3(k)=((RG/cpdet(Tplay(k)))
[1723]240     .       *((zsnet(k+1,nsza0(2),nlat0)-zsnet(k,nsza0(2),nlat0))
241     .         *PFRAC))
242     .      /((presrh(k)-presrh(k+1))*1.e5)
[1726]243       Qrh4(k)=((RG/cpdet(Tplay(k)))
[1723]244     .       *((zsnet(k+1,nsza0(2)-1,nlat0)-zsnet(k,nsza0(2)-1,nlat0))
245     .        *PFRAC))
246     .      /((presrh(k)-presrh(k+1))*1.e5)
247      ENDDO
[1726]248! Interapolation of PHEAT over GCM/MESOSCALE play levels
[1723]249      do j=1,klev
250        nl0 = nlrh-1
251        do i=nlrh,2,-1
[1726]252           if (exp(logplayrh(i)).ge.PPA(j)) then
[1723]253                nl0 = i-1
254           endif
255        enddo
256c        factflux = (log10(max(PPB(j),presrh(1)))-log10(presrh(nl0+1)))
257c     .            /(log10(presrh(nl0))-log10(presrh(nl0+1)))
[1726]258        factflux = (log10(max(PPA(j),exp(logplayrh(1))))
259     .                         -log10(exp(logplayrh(nl0+1))))
260     .     /(log10(exp(logplayrh(nl0)))-log10(exp(logplayrh(nl0+1))))
[1723]261        PHEATPPA(j)=factlat*(
[1726]262     .      factflux   *  factsza(2)  *Qrh3(nl0)
263     . +   factflux   *(1.-factsza(2))*Qrh4(nl0)
264     . + (1.-factflux)*  factsza(2)   *Qrh3(nl0+1)
265     . + (1.-factflux)*(1.-factsza(2))*Qrh4(nl0+1))
[1723]266     .            + (1.-factlat)*(
[1726]267     .      factflux   *  factsza(1)  *Qrh1(nl0)
268     . +   factflux   *(1.-factsza(1))*Qrh2(nl0)
269     . + (1.-factflux)*  factsza(1)   *Qrh1(nl0+1)
270     . + (1.-factflux)*(1.-factsza(1))*Qrh2(nl0+1) )
[1723]271        PHEAT(j)=PHEATPPA(j)
272      ENDDO
273
274
275#else     
[1591]276c Heating rates
277c -------------
278c On utilise le gradient du flux pour calculer le taux de chauffage:
279c   heat(K/s) = d(fluxnet)  (W/m2)
280c              *g           (m/s2)
281c              /(-dp)  (epaisseur couche, en Pa=kg/m/s2)
282c              /cp  (J/kg/K)
283
284      do j=1,klev
285! ADAPTATION GCM POUR CP(T)
[1726]286        PHEAT(j) = (ZFSNET(j+1)-ZFSNET(j))
[1591]287     .            *RG/cpdet(pt(j)) / ((PPB(j)-PPB(j+1))*1.e5)
288c-----TEST-------
289c tayloring the solar flux...
[2048]290c        if ((PPB(j).gt.0.04).and.(PPB(j).le.0.1)) then
291c         PHEAT(j) = PHEAT(j)*1.5
292c        endif
293c        if ((PPB(j).gt.0.1).and.(PPB(j).le.0.5)) then
294c         PHEAT(j) = PHEAT(j)*2.
295c        endif
[2464]296c BASE:
297c        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
298c          PHEAT(j) = PHEAT(j)*3
299c        endif
300c AFTER Tayloring TdeepD
301         if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
302           PHEAT(j) = PHEAT(j)*3.5
303         endif
304         if ((PPB(j).gt.10.).and.(PPB(j).le.50.)) then
305           PHEAT(j) = PHEAT(j)*1.5
306         endif
307c Options:
308c        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
309c        if ((PPB(j).gt.2.0).and.(PPB(j).le.10.)) then
[2048]310c        if ((PPB(j).gt.1.4).and.(PPB(j).le.100.)) then
[2464]311c          PHEAT(j) = PHEAT(j)*3.5
312c          PHEAT(j) = PHEAT(j)*3
313c          PHEAT(j) = PHEAT(j)*2.5
[2048]314c        endif
[2464]315c        if ((PPB(j).gt.10.).and.(PPB(j).le.35.)) then
316c        if ((PPB(j).gt.10.).and.(PPB(j).le.50.)) then
317c          PHEAT(j) = PHEAT(j)*2
318c          PHEAT(j) = PHEAT(j)*1.5
319c          PHEAT(j) = PHEAT(j)*1.3
320c        endif
321c        if ((PPB(j).gt.35.).and.(PPB(j).le.120.)) then
322c          PHEAT(j) = PHEAT(j)*2
323c          PHEAT(j) = PHEAT(j)*1.5
324c        endif
[1591]325c----------------
[1726]326      enddo
[1723]327#endif
328     
[1591]329
330      end
331
Note: See TracBrowser for help on using the repository browser.