source: trunk/LMDZ.VENUS/libf/phyvenus/radlwsw.F @ 1198

Last change on this file since 1198 was 973, checked in by slebonnois, 12 years ago

EM+SL: bug corrections in Venus physics
File size: 9.0 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $
3!
4      SUBROUTINE radlwsw(dist, rmu0, fract,
5     .                  paprs, pplay,tsol, t,
6     .                  heat,cool,radsol,
7     .                  topsw,toplw,solsw,sollw,
8     .                  sollwdown,
9     .                  lwnet, swnet)
10c     
11c======================================================================
12c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
13c Objet: interface entre le modele et les rayonnements
14c Arguments:
15c dist-----input-R- distance astronomique terre-soleil
16c rmu0-----input-R- cosinus de l'angle zenithal
17c fract----input-R- duree d'ensoleillement normalisee
18c solaire--input-R- constante solaire (W/m**2) (dans clesphys.h)
19c paprs----input-R- pression a inter-couche (Pa)
20c pplay----input-R- pression au milieu de couche (Pa)
21c tsol-----input-R- temperature du sol (en K)
22c t--------input-R- temperature (K)
23c heat-----output-R- echauffement atmospherique (visible) (K/jour)
24c cool-----output-R- refroidissement dans l'IR (K/jour)
25c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas)
26c topsw----output-R- flux solaire net au sommet de l'atm. (+ vers le bas)
27c toplw----output-R- ray. IR net au sommet de l'atmosphere (+ vers le haut)
28c solsw----output-R- flux solaire net a la surface (+ vers le bas)
29c sollw----output-R- ray. IR net a la surface (+ vers le bas)
30c sollwdown-output-R- ray. IR descendant a la surface (+ vers le bas)
31c lwnet____output-R- flux IR net (+ vers le haut)
32c swnet____output-R- flux solaire net (+ vers le bas)
33c
34     
35c MODIFS pour multimatrices ksi SPECIFIQUE VENUS
36c   S. Lebonnois    20/12/2006
37c   corrections     13/07/2007
38
39c======================================================================
40      use dimphy
41      USE comgeomphy
42      use write_field_phy
43      IMPLICIT none
44#include "dimensions.h"
45#include "YOMCST.h"
46#include "clesphys.h"
47#include "comcstVE.h"
48c
49      real rmu0(klon), fract(klon), dist
50c
51      real paprs(klon,klev+1), pplay(klon,klev)
52      real tsol(klon)
53      real t(klon,klev)
54      real heat(klon,klev), cool(klon,klev)
55      real radsol(klon), topsw(klon), toplw(klon)
56      real solsw(klon), sollw(klon)
57      real sollwdown(klon)
58      REAL swnet(klon,klev+1),lwnet(klon,klev+1)
59c
60      INTEGER k, kk, i, j, band
61c
62      REAL   PPB(klev+1)
63c
64      REAL   zfract, zrmu0
65c
66      REAL   zheat(klev), zcool(klev)
67      real  temp(klev)
68      REAL   ZFSNET(klev+1),ZFLNET(klev+1)
69      REAL   ztopsw, ztoplw
70      REAL   zsolsw, zsollw
71cIM BEG
72      REAL   zsollwdown
73cIM END
74      real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
75      real,save,allocatable :: ztop(:) ! in km
76
77      real    psi(0:klev+1,0:klev+1)
78      real    deltapsi(0:klev+1,0:klev+1)
79      real    latdeg
80      real    pt0(0:klev+1)
81      real    bplck(0:klev+1,nnuve)    ! Planck luminances in table layers
82      real    y(0:klev,nnuve)          ! intermediaire Planck
83      real    zdblay(0:klev+1,nnuve)   ! gradient en temperature de planck
84      integer mat,mat0
85      real    factp,factz,ksi
86
87      logical firstcall
88      data    firstcall/.true./
89      save    firstcall
90     
91c-------------------------------------------
92c  Initialisations
93c-----------------
94
95      if (firstcall) then
96
97c ---------- ksive --------------
98        allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
99        call load_ksi(ksive)
100c ---------- ztop --------------
101        allocate(ztop(klon))
102        DO i = 1, klon
103             ztop(i) = 70.
104        ENDDO !i
105c ztop: d'apres fit à figure 16 du papier Zavosa et al (tmp) traitant des
106c       donnees Venera
107c       DO i = 1, klon
108c         latdeg = abs(rlatd(i))
109c         if (latdeg.lt.15) then
110c            ztop(i) = 70.
111c         elseif (latdeg.lt.50) then
112c            ztop(i) = 63.95+6*cos((latdeg-15)*RPI/2./50.)
113c         else
114c            ztop(i) = min(63.95+6*cos((latdeg-15)*RPI/2./50.),
115c    .                     63.95-5.9*sin((latdeg-60)*RPI/2/30))
116c         endif
117c       print*,'lat(',i,')=',latdeg,'  ztop=',ztop(i)
118c       ENDDO !i
119c ---------- ztop --------------
120
121      endif ! firstcall
122c-------------------------------------------
123
124      DO k = 1, klev
125       DO i = 1, klon
126         heat(i,k)=0.
127         cool(i,k)=0.
128       ENDDO
129      ENDDO
130
131c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
132      DO j = 1, klon
133
134c======================================================================
135c  Initialisations
136c ---------------
137
138       DO k = 1, klev
139        zheat(k) = 0.0
140        zcool(k) = 0.0
141       ENDDO
142       DO k = 1, klev+1
143        ZFLNET(k) = 0.0
144        ZFSNET(k) = 0.0
145       ENDDO
146       ztopsw = 0.0
147       ztoplw = 0.0
148       zsolsw = 0.0
149       zsollw = 0.0
150       zsollwdown = 0.0
151     
152         zfract = fract(j)
153         zrmu0 = rmu0(j)
154 
155       DO k = 1, klev+1
156         PPB(k) = paprs(j,k)/1.e5
157       ENDDO
158
159       pt0(0)  = tsol(j)
160       DO k = 1, klev
161         pt0(k) = t(j,k)
162       ENDDO
163       pt0(klev+1) = 0.
164 
165       DO k = 0,klev+1
166       DO i = 0,klev+1
167        psi(i,k) = 0.   ! positif quand nrj de i->k
168        deltapsi(i,k) = 0.
169       ENDDO
170       ENDDO
171       
172c======================================================================
173c Getting psi and deltapsi
174c ------------------------
175
176c Planck function
177c ---------------
178      do band=1,nnuve
179        do k=0,klev
180c B(T,l) = al/(exp(bl/T)-1)
181         y(k,band) = exp(bl(band)/pt0(k))-1.
182         bplck(k,band) = al(band)/(y(k,band))
183         zdblay(k,band)= al(band)*bl(band)*exp(bl(band)/pt0(k))/
184     .                  ((pt0(k)*pt0(k))*(y(k,band)*y(k,band)))
185        enddo
186        bplck(klev+1,band) = 0.0
187        zdblay(klev+1,band)= 0.0
188      enddo
189
190c finding the right matrixes
191c --------------------------
192       mat0 = 0
193       do mat=1,nbmat-nbztopve
194         if (  (psurfve(mat).ge.paprs(j,1))
195     .    .and.(psurfve(mat+nbztopve).lt.paprs(j,1))
196     .    .and.(ztopve(mat).lt.ztop(j))
197     .    .and.(ztopve(mat+1).ge.ztop(j)) ) then
198              mat0  = mat
199c             print*,'ig=',j,'  mat0=',mat
200              factp = (paprs(j,1)           -psurfve(mat))
201     .               /(psurfve(mat+nbztopve)-psurfve(mat))
202              factz = (ztop(j)      -ztopve(mat))
203     .               /(ztopve(mat+1)-ztopve(mat))
204              exit
205         endif
206       enddo
207       if (mat0.eq.0) then
208         write(*,*) 'Finding the right matrix in radlwsw'
209         print*,'Probleme pour interpolation au point ig=',j
210         print*,'psurf = ',paprs(j,1),' ztop = ',ztop(j)
211         stop
212       endif
213
214c interpolation of ksi and computation of psi,deltapsi
215c ----------------------------------------------------
216       do band=1,nnuve
217        do k=0,klev+1
218         do i=0,klev+1
219          ksi = ksive(i,k,band,mat0)*(1-factz)*(1-factp)
220     .         +ksive(i,k,band,mat0+1)*factz  *(1-factp)
221     .         +ksive(i,k,band,mat0+nbztopve)*(1-factz)*factp
222     .         +ksive(i,k,band,mat0+nbztopve+1)*factz  *factp
223          psi(i,k) = psi(i,k) +
224     .               ksi*(bplck(i,band)-bplck(k,band))
225          deltapsi(i,k) = deltapsi(i,k) + ksi*zdblay(i,band)
226         enddo
227        enddo
228       enddo
229
230c======================================================================
231c LW call
232c---------
233      temp(1:klev)=t(j,1:klev)
234      CALL LW_venus_ve(
235     .        PPB,temp,psi,deltapsi,
236     .        zcool,
237     .        ztoplw,zsollw,
238     .        zsollwdown,ZFLNET)
239
240c---------
241c SW call
242c---------
243      CALL SW_venus_dc(zrmu0, zfract,
244     S        PPB,temp,
245     S        zheat,
246     S        ztopsw,zsolsw,ZFSNET)
247     
248c======================================================================
249         radsol(j) = zsolsw - zsollw  ! + vers bas
250         topsw(j) = ztopsw            ! + vers bas
251         toplw(j) = ztoplw            ! + vers haut
252         solsw(j) = zsolsw            ! + vers bas
253         sollw(j) = -zsollw           ! + vers bas
254         sollwdown(j) = zsollwdown    ! + vers bas
255
256         DO k = 1, klev+1
257         lwnet  (j,k)   = ZFLNET(k)
258         swnet  (j,k)   = ZFSNET(k)
259         ENDDO
260
261      DO k = 1, klev
262         heat (j,k) = zheat(k)
263         cool (j,k) = zcool(k)
264      ENDDO
265c
266      ENDDO ! of DO j = 1, klon
267c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
268! for tests: write output fields...
269!      call writefield_phy('radlwsw_heat',heat,klev)
270!      call writefield_phy('radlwsw_cool',cool,klev)
271!      call writefield_phy('radlwsw_radsol',radsol,1)
272!      call writefield_phy('radlwsw_topsw',topsw,1)
273!      call writefield_phy('radlwsw_toplw',toplw,1)
274!      call writefield_phy('radlwsw_solsw',solsw,1)
275!      call writefield_phy('radlwsw_sollw',sollw,1)
276!      call writefield_phy('radlwsw_sollwdown',sollwdown,1)
277!      call writefield_phy('radlwsw_swnet',swnet,klev+1)
278!      call writefield_phy('radlwsw_lwnet',lwnet,klev+1)
279
280c tests
281
282c     j = klon/2
283c     j = 1
284c     print*,'mu0=',rmu0(j)
285c     print*,'   net flux vis   HEAT(K/day)'
286c     do k=1,klev
287c     print*,k,ZFSNET(k),heat(j,k)*8.56548e-3
288c     enddo
289c     print*,'   net flux IR    COOL(K/day)'
290c     do k=1,klev
291c     print*,k,ZFLNET(k),cool(j,k)*8.56548e-3
292c     enddo
293
294      firstcall = .false.
295      RETURN
296      END
297
Note: See TracBrowser for help on using the repository browser.