! ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $ ! SUBROUTINE radlwsw(dist, rmu0, fract, . paprs, pplay,tsol, t, . heat,cool,radsol, . topsw,toplw,solsw,sollw, . sollwdown, . lwnet, swnet) c c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719 c Objet: interface entre le modele et les rayonnements c Arguments: c dist-----input-R- distance astronomique terre-soleil c rmu0-----input-R- cosinus de l'angle zenithal c fract----input-R- duree d'ensoleillement normalisee c solaire--input-R- constante solaire (W/m**2) (dans clesphys.h) c paprs----input-R- pression a inter-couche (Pa) c pplay----input-R- pression au milieu de couche (Pa) c tsol-----input-R- temperature du sol (en K) c t--------input-R- temperature (K) c heat-----output-R- echauffement atmospherique (visible) (K/jour) c cool-----output-R- refroidissement dans l'IR (K/jour) c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas) c topsw----output-R- flux solaire net au sommet de l'atm. (+ vers le bas) c toplw----output-R- ray. IR net au sommet de l'atmosphere (+ vers le haut) c solsw----output-R- flux solaire net a la surface (+ vers le bas) c sollw----output-R- ray. IR net a la surface (+ vers le bas) c sollwdown-output-R- ray. IR descendant a la surface (+ vers le bas) c lwnet____output-R- flux IR net (+ vers le haut) c swnet____output-R- flux solaire net (+ vers le bas) c c MODIFS pour multimatrices ksi SPECIFIQUE VENUS c S. Lebonnois 20/12/2006 c corrections 13/07/2007 c====================================================================== use dimphy USE comgeomphy use write_field_phy IMPLICIT none #include "dimensions.h" #include "YOMCST.h" #include "clesphys.h" #include "comcstVE.h" c real rmu0(klon), fract(klon), dist c real paprs(klon,klev+1), pplay(klon,klev) real tsol(klon) real t(klon,klev) real heat(klon,klev), cool(klon,klev) real radsol(klon), topsw(klon), toplw(klon) real solsw(klon), sollw(klon) real sollwdown(klon) REAL swnet(klon,klev+1),lwnet(klon,klev+1) c INTEGER k, kk, i, j, band c REAL PPB(klev+1) c REAL zfract, zrmu0 c REAL zheat(klev), zcool(klev) real temp(klev) REAL ZFSNET(klev+1),ZFLNET(klev+1) REAL ztopsw, ztoplw REAL zsolsw, zsollw cIM BEG REAL zsollwdown cIM END real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file real,save,allocatable :: ztop(:) ! in km real psi(0:klev+1,0:klev+1) real deltapsi(0:klev+1,0:klev+1) real latdeg real pt0(0:klev+1) real bplck(0:klev+1,nnuve) ! Planck luminances in table layers real y(0:klev,nnuve) ! intermediaire Planck real zdblay(0:klev+1,nnuve) ! gradient en temperature de planck integer mat,mat0 real factp,factz,ksi logical firstcall data firstcall/.true./ save firstcall c------------------------------------------- c Initialisations c----------------- if (firstcall) then c ---------- ksive -------------- allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat)) call load_ksi(ksive) c ---------- ztop -------------- allocate(ztop(klon)) DO i = 1, klon ztop(i) = 70. ENDDO !i c ztop: d'apres fit à figure 16 du papier Zavosa et al (tmp) traitant des c donnees Venera c DO i = 1, klon c latdeg = abs(rlatd(i)) c if (latdeg.lt.15) then c ztop(i) = 70. c elseif (latdeg.lt.50) then c ztop(i) = 63.95+6*cos((latdeg-15)*RPI/2./50.) c else c ztop(i) = min(63.95+6*cos((latdeg-15)*RPI/2./50.), c . 63.95-5.9*sin((latdeg-60)*RPI/2/30)) c endif c print*,'lat(',i,')=',latdeg,' ztop=',ztop(i) c ENDDO !i c ---------- ztop -------------- endif ! firstcall c------------------------------------------- DO k = 1, klev DO i = 1, klon heat(i,k)=0. cool(i,k)=0. ENDDO ENDDO c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++ DO j = 1, klon c====================================================================== c Initialisations c --------------- DO k = 1, klev zheat(k) = 0.0 zcool(k) = 0.0 ENDDO DO k = 1, klev+1 ZFLNET(k) = 0.0 ZFSNET(k) = 0.0 ENDDO ztopsw = 0.0 ztoplw = 0.0 zsolsw = 0.0 zsollw = 0.0 zsollwdown = 0.0 zfract = fract(j) zrmu0 = rmu0(j) DO k = 1, klev+1 PPB(k) = paprs(j,k)/1.e5 ENDDO pt0(0) = tsol(j) DO k = 1, klev pt0(k) = t(j,k) ENDDO pt0(klev+1) = 0. DO k = 0,klev+1 DO i = 0,klev+1 psi(i,k) = 0. ! positif quand nrj de i->k deltapsi(i,k) = 0. ENDDO ENDDO c====================================================================== c Getting psi and deltapsi c ------------------------ c Planck function c --------------- do band=1,nnuve do k=0,klev c B(T,l) = al/(exp(bl/T)-1) y(k,band) = exp(bl(band)/pt0(k))-1. bplck(k,band) = al(band)/(y(k,band)) zdblay(k,band)= al(band)*bl(band)*exp(bl(band)/pt0(k))/ . ((pt0(k)*pt0(k))*(y(k,band)*y(k,band))) enddo bplck(klev+1,band) = 0.0 zdblay(klev+1,band)= 0.0 enddo c finding the right matrixes c -------------------------- mat0 = 0 do mat=1,nbmat-nbztopve if ( (psurfve(mat).ge.paprs(j,1)) . .and.(psurfve(mat+nbztopve).lt.paprs(j,1)) . .and.(ztopve(mat).lt.ztop(j)) . .and.(ztopve(mat+1).ge.ztop(j)) ) then mat0 = mat c print*,'ig=',j,' mat0=',mat factp = (paprs(j,1) -psurfve(mat)) . /(psurfve(mat+nbztopve)-psurfve(mat)) factz = (ztop(j) -ztopve(mat)) . /(ztopve(mat+1)-ztopve(mat)) exit endif enddo if (mat0.eq.0) then write(*,*) 'Finding the right matrix in radlwsw' print*,'Probleme pour interpolation au point ig=',j print*,'psurf = ',paprs(j,1),' ztop = ',ztop(j) stop endif c interpolation of ksi and computation of psi,deltapsi c ---------------------------------------------------- do band=1,nnuve do k=0,klev+1 do i=0,klev+1 ksi = ksive(i,k,band,mat0)*(1-factz)*(1-factp) . +ksive(i,k,band,mat0+1)*factz *(1-factp) . +ksive(i,k,band,mat0+nbztopve)*(1-factz)*factp . +ksive(i,k,band,mat0+nbztopve+1)*factz *factp psi(i,k) = psi(i,k) + . ksi*(bplck(i,band)-bplck(k,band)) deltapsi(i,k) = deltapsi(i,k) + ksi*zdblay(i,band) enddo enddo enddo c====================================================================== c LW call c--------- temp(1:klev)=t(j,1:klev) CALL LW_venus_ve( . PPB,temp,psi,deltapsi, . zcool, . ztoplw,zsollw, . zsollwdown,ZFLNET) c--------- c SW call c--------- CALL SW_venus_dc(zrmu0, zfract, S PPB,temp, S zheat, S ztopsw,zsolsw,ZFSNET) c====================================================================== radsol(j) = zsolsw - zsollw ! + vers bas topsw(j) = ztopsw ! + vers bas toplw(j) = ztoplw ! + vers haut solsw(j) = zsolsw ! + vers bas sollw(j) = -zsollw ! + vers bas sollwdown(j) = zsollwdown ! + vers bas DO k = 1, klev+1 lwnet (j,k) = ZFLNET(k) swnet (j,k) = ZFSNET(k) ENDDO DO k = 1, klev heat (j,k) = zheat(k) cool (j,k) = zcool(k) ENDDO c ENDDO ! of DO j = 1, klon c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++ ! for tests: write output fields... ! call writefield_phy('radlwsw_heat',heat,klev) ! call writefield_phy('radlwsw_cool',cool,klev) ! call writefield_phy('radlwsw_radsol',radsol,1) ! call writefield_phy('radlwsw_topsw',topsw,1) ! call writefield_phy('radlwsw_toplw',toplw,1) ! call writefield_phy('radlwsw_solsw',solsw,1) ! call writefield_phy('radlwsw_sollw',sollw,1) ! call writefield_phy('radlwsw_sollwdown',sollwdown,1) ! call writefield_phy('radlwsw_swnet',swnet,klev+1) ! call writefield_phy('radlwsw_lwnet',lwnet,klev+1) c tests c j = klon/2 c j = 1 c print*,'mu0=',rmu0(j) c print*,' net flux vis HEAT(K/day)' c do k=1,klev c print*,k,ZFSNET(k),heat(j,k)*8.56548e-3 c enddo c print*,' net flux IR COOL(K/day)' c do k=1,klev c print*,k,ZFLNET(k),cool(j,k)*8.56548e-3 c enddo firstcall = .false. RETURN END