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

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

EM+SL: bug corrections in Venus physics
File size: 9.0 KB
RevLine 
[3]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======================================================================
[101]40      use dimphy
41      USE comgeomphy
[973]42      use write_field_phy
[101]43      IMPLICIT none
[3]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)
[892]58      REAL swnet(klon,klev+1),lwnet(klon,klev+1)
[3]59c
[953]60      INTEGER k, kk, i, j, band
[3]61c
[892]62      REAL   PPB(klev+1)
[3]63c
64      REAL   zfract, zrmu0
65c
[892]66      REAL   zheat(klev), zcool(klev)
[973]67      real  temp(klev)
[892]68      REAL   ZFSNET(klev+1),ZFLNET(klev+1)
[3]69      REAL   ztopsw, ztoplw
70      REAL   zsolsw, zsollw
71cIM BEG
72      REAL   zsollwdown
73cIM END
[101]74      real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
[953]75      real,save,allocatable :: ztop(:) ! in km
76
[892]77      real    psi(0:klev+1,0:klev+1)
78      real    deltapsi(0:klev+1,0:klev+1)
[101]79      real    latdeg
[953]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
[3]86
87      logical firstcall
88      data    firstcall/.true./
89      save    firstcall
90     
91c-------------------------------------------
92c  Initialisations
93c-----------------
94
95      if (firstcall) then
[101]96
97c ---------- ksive --------------
[892]98        allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
[3]99        call load_ksi(ksive)
100c ---------- ztop --------------
[101]101        allocate(ztop(klon))
[3]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
[953]122c-------------------------------------------
[3]123
124      DO k = 1, klev
[953]125       DO i = 1, klon
[3]126         heat(i,k)=0.
127         cool(i,k)=0.
[953]128       ENDDO
[3]129      ENDDO
[953]130
[3]131c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
[973]132      DO j = 1, klon
[953]133
134c======================================================================
135c  Initialisations
136c ---------------
137
[3]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 
[953]155       DO k = 1, klev+1
[3]156         PPB(k) = paprs(j,k)/1.e5
[953]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.
[3]164 
[953]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
[3]171       
172c======================================================================
[953]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======================================================================
[3]231c LW call
232c---------
[973]233      temp(1:klev)=t(j,1:klev)
[3]234      CALL LW_venus_ve(
[973]235     .        PPB,temp,psi,deltapsi,
[3]236     .        zcool,
237     .        ztoplw,zsollw,
238     .        zsollwdown,ZFLNET)
239
240c---------
241c SW call
242c---------
243      CALL SW_venus_dc(zrmu0, zfract,
[973]244     S        PPB,temp,
[3]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
[892]256         DO k = 1, klev+1
[3]257         lwnet  (j,k)   = ZFLNET(k)
258         swnet  (j,k)   = ZFSNET(k)
259         ENDDO
260
[892]261      DO k = 1, klev
[3]262         heat (j,k) = zheat(k)
263         cool (j,k) = zcool(k)
264      ENDDO
265c
[973]266      ENDDO ! of DO j = 1, klon
[3]267c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
[973]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)
[3]279
280c tests
281
282c     j = klon/2
283c     j = 1
284c     print*,'mu0=',rmu0(j)
285c     print*,'   net flux vis   HEAT(K/day)'
[892]286c     do k=1,klev
[3]287c     print*,k,ZFSNET(k),heat(j,k)*8.56548e-3
288c     enddo
289c     print*,'   net flux IR    COOL(K/day)'
[892]290c     do k=1,klev
[3]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.