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

!
! $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
      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, nb_gr
c
      REAL   PPB(klev+1)
c
      REAL   zfract, zrmu0
c
      REAL   zheat(klev), zcool(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    psimap(0:klev+1,0:klev+1,klon)
      real    deltapsimap(0:klev+1,0:klev+1,klon)
      real    psi(0:klev+1,0:klev+1)
      real    deltapsi(0:klev+1,0:klev+1)
      real    latdeg
      real    pt0(klon,0:klev+1)
      real,save,allocatable :: ztop(:) ! in km

      logical firstcall
      data    firstcall/.true./
      save    firstcall
      
c-------------------------------------------
      nb_gr = klon
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

      DO i = 1, klon
          pt0(i,0)  = tsol(i)
          DO k = 1, klev
            pt0(i,k) = t(i,k)
          ENDDO
          pt0(i,klev+1) = 0.
      ENDDO !i

      call load_psi(paprs(:,1),ztop,ksive,pt0,psimap,deltapsimap)

      DO k = 1, klev
      DO i = 1, klon
         heat(i,k)=0.
         cool(i,k)=0.
      ENDDO
      ENDDO
c
c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
      DO 99999 j = 1, nb_gr
 
       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
 
      DO k = 0,klev+1
      DO i = 0,klev+1
        psi(i,k) = psimap(i,k,j) 
        deltapsi(i,k) = deltapsimap(i,k,j) 
      ENDDO
      ENDDO
       
c======================================================================
c LW call
c---------
      CALL LW_venus_ve(
     .        PPB,t(j,:),psi,deltapsi,
     .        zcool,
     .        ztoplw,zsollw,
     .        zsollwdown,ZFLNET)

c---------
c SW call
c---------
      CALL SW_venus_dc(zrmu0, zfract,
     S        PPB,t(j,:), 
     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
99999 CONTINUE
c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++

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