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

Last change on this file since 963 was 953, checked in by slebonnois, 12 years ago

SL: optimisation pour le parallèle suite à tests Venus / petite correction appels routines secondaires dans Venus et Titan
File size: 8.4 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      IMPLICIT none
43#include "dimensions.h"
44#include "YOMCST.h"
45#include "clesphys.h"
46#include "comcstVE.h"
47c
48      real rmu0(klon), fract(klon), dist
49c
50      real paprs(klon,klev+1), pplay(klon,klev)
51      real tsol(klon)
52      real t(klon,klev)
53      real heat(klon,klev), cool(klon,klev)
54      real radsol(klon), topsw(klon), toplw(klon)
55      real solsw(klon), sollw(klon)
56      real sollwdown(klon)
57      REAL swnet(klon,klev+1),lwnet(klon,klev+1)
58c
59      INTEGER k, kk, i, j, band
60c
61      REAL   PPB(klev+1)
62c
63      REAL   zfract, zrmu0
64c
65      REAL   zheat(klev), zcool(klev)
66      REAL   ZFSNET(klev+1),ZFLNET(klev+1)
67      REAL   ztopsw, ztoplw
68      REAL   zsolsw, zsollw
69cIM BEG
70      REAL   zsollwdown
71cIM END
72      real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
73      real,save,allocatable :: ztop(:) ! in km
74
75      real    psi(0:klev+1,0:klev+1)
76      real    deltapsi(0:klev+1,0:klev+1)
77      real    latdeg
78      real    pt0(0:klev+1)
79      real    bplck(0:klev+1,nnuve)    ! Planck luminances in table layers
80      real    y(0:klev,nnuve)          ! intermediaire Planck
81      real    zdblay(0:klev+1,nnuve)   ! gradient en temperature de planck
82      integer mat,mat0
83      real    factp,factz,ksi
84
85      logical firstcall
86      data    firstcall/.true./
87      save    firstcall
88     
89c-------------------------------------------
90c  Initialisations
91c-----------------
92
93      if (firstcall) then
94
95c ---------- ksive --------------
96        allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
97        call load_ksi(ksive)
98c ---------- ztop --------------
99        allocate(ztop(klon))
100        DO i = 1, klon
101             ztop(i) = 70.
102        ENDDO !i
103c ztop: d'apres fit à figure 16 du papier Zavosa et al (tmp) traitant des
104c       donnees Venera
105c       DO i = 1, klon
106c         latdeg = abs(rlatd(i))
107c         if (latdeg.lt.15) then
108c            ztop(i) = 70.
109c         elseif (latdeg.lt.50) then
110c            ztop(i) = 63.95+6*cos((latdeg-15)*RPI/2./50.)
111c         else
112c            ztop(i) = min(63.95+6*cos((latdeg-15)*RPI/2./50.),
113c    .                     63.95-5.9*sin((latdeg-60)*RPI/2/30))
114c         endif
115c       print*,'lat(',i,')=',latdeg,'  ztop=',ztop(i)
116c       ENDDO !i
117c ---------- ztop --------------
118
119      endif ! firstcall
120c-------------------------------------------
121
122      DO k = 1, klev
123       DO i = 1, klon
124         heat(i,k)=0.
125         cool(i,k)=0.
126       ENDDO
127      ENDDO
128
129c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
130      DO 99999 j = 1, klon
131
132c======================================================================
133c  Initialisations
134c ---------------
135
136       DO k = 1, klev
137        zheat(k) = 0.0
138        zcool(k) = 0.0
139       ENDDO
140       DO k = 1, klev+1
141        ZFLNET(k) = 0.0
142        ZFSNET(k) = 0.0
143       ENDDO
144       ztopsw = 0.0
145       ztoplw = 0.0
146       zsolsw = 0.0
147       zsollw = 0.0
148       zsollwdown = 0.0
149     
150         zfract = fract(j)
151         zrmu0 = rmu0(j)
152 
153       DO k = 1, klev+1
154         PPB(k) = paprs(j,k)/1.e5
155       ENDDO
156
157       pt0(0)  = tsol(j)
158       DO k = 1, klev
159         pt0(k) = t(j,k)
160       ENDDO
161       pt0(klev+1) = 0.
162 
163       DO k = 0,klev+1
164       DO i = 0,klev+1
165        psi(i,k) = 0.   ! positif quand nrj de i->k
166        deltapsi(i,k) = 0.
167       ENDDO
168       ENDDO
169       
170c======================================================================
171c Getting psi and deltapsi
172c ------------------------
173
174c Planck function
175c ---------------
176      do band=1,nnuve
177        do k=0,klev
178c B(T,l) = al/(exp(bl/T)-1)
179         y(k,band) = exp(bl(band)/pt0(k))-1.
180         bplck(k,band) = al(band)/(y(k,band))
181         zdblay(k,band)= al(band)*bl(band)*exp(bl(band)/pt0(k))/
182     .                  ((pt0(k)*pt0(k))*(y(k,band)*y(k,band)))
183        enddo
184        bplck(klev+1,band) = 0.0
185        zdblay(klev+1,band)= 0.0
186      enddo
187
188c finding the right matrixes
189c --------------------------
190       mat0 = 0
191       do mat=1,nbmat-nbztopve
192         if (  (psurfve(mat).ge.paprs(j,1))
193     .    .and.(psurfve(mat+nbztopve).lt.paprs(j,1))
194     .    .and.(ztopve(mat).lt.ztop(j))
195     .    .and.(ztopve(mat+1).ge.ztop(j)) ) then
196              mat0  = mat
197c             print*,'ig=',j,'  mat0=',mat
198              factp = (paprs(j,1)           -psurfve(mat))
199     .               /(psurfve(mat+nbztopve)-psurfve(mat))
200              factz = (ztop(j)      -ztopve(mat))
201     .               /(ztopve(mat+1)-ztopve(mat))
202              exit
203         endif
204       enddo
205       if (mat0.eq.0) then
206         write(*,*) 'Finding the right matrix in radlwsw'
207         print*,'Probleme pour interpolation au point ig=',j
208         print*,'psurf = ',paprs(j,1),' ztop = ',ztop(j)
209         stop
210       endif
211
212c interpolation of ksi and computation of psi,deltapsi
213c ----------------------------------------------------
214       do band=1,nnuve
215        do k=0,klev+1
216         do i=0,klev+1
217          ksi = ksive(i,k,band,mat0)*(1-factz)*(1-factp)
218     .         +ksive(i,k,band,mat0+1)*factz  *(1-factp)
219     .         +ksive(i,k,band,mat0+nbztopve)*(1-factz)*factp
220     .         +ksive(i,k,band,mat0+nbztopve+1)*factz  *factp
221          psi(i,k) = psi(i,k) +
222     .               ksi*(bplck(i,band)-bplck(k,band))
223          deltapsi(i,k) = deltapsi(i,k) + ksi*zdblay(i,band)
224         enddo
225        enddo
226       enddo
227
228c======================================================================
229c LW call
230c---------
231      CALL LW_venus_ve(
232     .        PPB,t(j,:),psi,deltapsi,
233     .        zcool,
234     .        ztoplw,zsollw,
235     .        zsollwdown,ZFLNET)
236
237c---------
238c SW call
239c---------
240      CALL SW_venus_dc(zrmu0, zfract,
241     S        PPB,t(j,:),
242     S        zheat,
243     S        ztopsw,zsolsw,ZFSNET)
244     
245c======================================================================
246         radsol(j) = zsolsw - zsollw  ! + vers bas
247         topsw(j) = ztopsw            ! + vers bas
248         toplw(j) = ztoplw            ! + vers haut
249         solsw(j) = zsolsw            ! + vers bas
250         sollw(j) = -zsollw           ! + vers bas
251         sollwdown(j) = zsollwdown    ! + vers bas
252
253         DO k = 1, klev+1
254         lwnet  (j,k)   = ZFLNET(k)
255         swnet  (j,k)   = ZFSNET(k)
256         ENDDO
257
258      DO k = 1, klev
259         heat (j,k) = zheat(k)
260         cool (j,k) = zcool(k)
261      ENDDO
262c
26399999 CONTINUE
264c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
265
266c tests
267
268c     j = klon/2
269c     j = 1
270c     print*,'mu0=',rmu0(j)
271c     print*,'   net flux vis   HEAT(K/day)'
272c     do k=1,klev
273c     print*,k,ZFSNET(k),heat(j,k)*8.56548e-3
274c     enddo
275c     print*,'   net flux IR    COOL(K/day)'
276c     do k=1,klev
277c     print*,k,ZFLNET(k),cool(j,k)*8.56548e-3
278c     enddo
279
280      firstcall = .false.
281      RETURN
282      END
283
Note: See TracBrowser for help on using the repository browser.