source: LMDZ6/branches/Amaury_dev/libf/phylmd/simu_airs.F90 @ 5093

        
        module m_simu_airs

        USE print_control_mod, ONLY: prt_level,lunout
          
        implicit none

        REAL, PARAMETER :: tau_thresh = 0.05 ! seuil nuages detectables
        REAL, PARAMETER :: p_thresh = 445.   ! seuil nuages hauts
        REAL, PARAMETER :: em_min = 0.2      ! seuils nuages semi-transp
        REAL, PARAMETER :: em_max = 0.85
        REAL, PARAMETER :: undef = 999.

        contains

        REAL function search_tropopause(P,T,alt,N) result(P_tropo)
! this function searches for the tropopause pressure in [hPa].
! The search is based on ideology described in
! Reichler et al., Determining the tropopause height from gridded data,
! GRL, 30(20) 2042, doi:10.1029/2003GL018240, 2003

        INTEGER N,i,i_lev,first_point,exit_flag,i_dir
        REAL P(N),T(N),alt(N),slope(N)
        REAL P_min, P_max, slope_limit,slope_2km, &
   delta_alt_limit,tmp,delta_alt
        PARAMETER(P_min=75.0, P_max=470.0)   ! hPa
        PARAMETER(slope_limit=0.002)         ! 2 K/km converted to K/m
        PARAMETER(delta_alt_limit=2000.0)    ! 2000 meters

        do i=1,N-1
        slope(i)=-(T(i+1)-T(i))/(alt(i+1)-alt(i))
        END DO
        slope(N)=slope(N-1)

        if (P(1)>P(N)) then
        i_dir= 1
        i=1
        else
        i_dir=-1
        i=N
        end if

        first_point=0
        exit_flag=0
        do while(exit_flag==0)
        if (P(i)>=P_min.and.P(i)<=P_max) then
        if (first_point>0) then
        delta_alt=alt(i)-alt(first_point)
        if (delta_alt>=delta_alt_limit) then
        slope_2km=(T(first_point)-T(i))/delta_alt
        if (slope_2km<slope_limit) then
        exit_flag=1
        else
        first_point=0
        end if
        end if
        end if
        if (first_point==0.and.slope(i)<slope_limit) first_point=i
        end if
        i=i+i_dir
        if (i<=1.or.i>=N) exit_flag=1
        END DO

        if (first_point<=0) P_tropo=65.4321
        if (first_point==1) P_tropo=64.5432
        if (first_point>1) then
        tmp=(slope_limit-slope(first_point))/(slope(first_point+1)- &
   slope(first_point))*(P(first_point+1)-P(first_point))
        P_tropo=P(first_point)+tmp
        ! print*, 'P_tropo= ', tmp, P(first_point), P_tropo
        end if

! Ajout Marine
        if (P_tropo < 60. .or. P_tropo > 470.) then
        P_tropo = 999.
        endif

        return
        end function search_tropopause



        subroutine cloud_structure(len_cs, rneb_cs, temp_cs, &
   emis_cs, iwco_cs, &
   pres_cs, dz_cs, rhodz_cs, rad_cs, &
   cc_tot_cs, cc_hc_cs, cc_hist_cs, &
   cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
   pcld_hc_cs, tcld_hc_cs, &
   em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
   pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
   pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
   pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
   em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)


     
        INTEGER :: i, n, nss

        INTEGER, intent(in) :: len_cs
        REAL, DIMENSION(:), intent(in) :: rneb_cs, temp_cs
        REAL, DIMENSION(:), intent(in) :: emis_cs, iwco_cs, rad_cs
        REAL, DIMENSION(:), intent(in) :: pres_cs, dz_cs, rhodz_cs

        REAL, intent(out) :: cc_tot_cs, cc_hc_cs, cc_hist_cs, &
   cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
   pcld_hc_cs, tcld_hc_cs, em_hc_cs, iwp_hc_cs, &
   pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
   pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
   pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
   em_hist_cs, iwp_hist_cs, &
   deltaz_hc_cs, deltaz_hist_cs, rad_hist_cs

        REAL, DIMENSION(len_cs) :: rneb_ord
        REAL :: rneb_min

        REAL, DIMENSION(:), allocatable :: s, s_hc, s_hist, rneb_max
        REAL, DIMENSION(:), allocatable :: sCb, sThCi, sAnv
        REAL, DIMENSION(:), allocatable :: iwp_ss, pcld_ss, tcld_ss,&
   emis_ss
        REAL, DIMENSION(:), allocatable :: deltaz_ss, rad_ss

        CHARACTER (len = 50)      :: modname = 'simu_airs.cloud_structure'
        CHARACTER (len = 160)     :: abort_message
        

        write(lunout,*) 'dans cloud_structure'

        call ordonne(len_cs, rneb_cs, rneb_ord)
       

! Definition des sous_sections

        rneb_min = rneb_ord(1)
        nss = 1

        do i = 2, size(rneb_cs)
        if (rneb_ord(i) > rneb_min) then
        nss = nss+1
        rneb_min = rneb_ord(i)
        endif
        enddo

        allocate (s(nss))
        allocate (s_hc(nss))
        allocate (s_hist(nss))
        allocate (rneb_max(nss))
        allocate (emis_ss(nss))
        allocate (pcld_ss(nss))
        allocate (tcld_ss(nss))
        allocate (iwp_ss(nss))
        allocate (deltaz_ss(nss))
        allocate (rad_ss(nss))
        allocate (sCb(nss))
        allocate (sThCi(nss))
        allocate (sAnv(nss))

        rneb_min = rneb_ord(1)
        n = 1
        s(1) = rneb_ord(1)
        s_hc(1) = rneb_ord(1)
        s_hist(1) = rneb_ord(1)
        sCb(1) = rneb_ord(1)
        sThCi(1) = rneb_ord(1)
        sAnv(1) = rneb_ord(1)

        rneb_max(1) = rneb_ord(1)

        do i = 2, size(rneb_cs)
        if (rneb_ord(i) > rneb_min) then
        n = n+1
        s(n) = rneb_ord(i)-rneb_min
        s_hc(n) = rneb_ord(i)-rneb_min
        s_hist(n) = rneb_ord(i)-rneb_min
        sCb(n) = rneb_ord(i)-rneb_min
        sThCi(n) = rneb_ord(i)-rneb_min
        sAnv(n) = rneb_ord(i)-rneb_min

        rneb_max(n) = rneb_ord(i)
        rneb_min = rneb_ord(i)
        endif
        enddo

! Appel de sous_section

        do i = 1, nss
         call sous_section(len_cs, rneb_cs, temp_cs, &
    emis_cs, iwco_cs, &
    pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
    rneb_max(i),s(i),s_hc(i),s_hist(i), &
    sCb(i), sThCi(i), sAnv(i), &
    emis_ss(i), &
    pcld_ss(i), tcld_ss(i), iwp_ss(i), deltaz_ss(i), rad_ss(i))
        enddo

! Caracteristiques de la structure nuageuse

        cc_tot_cs = 0.
        cc_hc_cs = 0.
        cc_hist_cs = 0.

        cc_Cb_cs = 0.
        cc_ThCi_cs = 0.
        cc_Anv_cs = 0.

        em_hc_cs = 0.
        iwp_hc_cs = 0.
        deltaz_hc_cs = 0.

        em_hist_cs = 0.
        iwp_hist_cs = 0.
        deltaz_hist_cs = 0.
        rad_hist_cs = 0.

        pcld_hc_cs = 0.
        tcld_hc_cs = 0.

        pcld_Cb_cs = 0.
        tcld_Cb_cs = 0.
        em_Cb_cs = 0.

        pcld_ThCi_cs = 0.
        tcld_ThCi_cs = 0.
        em_ThCi_cs = 0.

        pcld_Anv_cs = 0.
        tcld_Anv_cs = 0.
        em_Anv_cs = 0.

         do i = 1, nss

        cc_tot_cs = cc_tot_cs + s(i)
        cc_hc_cs = cc_hc_cs + s_hc(i)
        cc_hist_cs = cc_hist_cs + s_hist(i)

        cc_Cb_cs = cc_Cb_cs + sCb(i)
        cc_ThCi_cs = cc_ThCi_cs + sThCi(i)
        cc_Anv_cs = cc_Anv_cs + sAnv(i)

        iwp_hc_cs = iwp_hc_cs + s_hc(i)*iwp_ss(i)
        em_hc_cs = em_hc_cs + s_hc(i)*emis_ss(i)
        deltaz_hc_cs = deltaz_hc_cs + s_hc(i)*deltaz_ss(i)

        iwp_hist_cs = iwp_hist_cs + s_hist(i)*iwp_ss(i)
        em_hist_cs = em_hist_cs + s_hist(i)*emis_ss(i)
        deltaz_hist_cs = deltaz_hist_cs + s_hist(i)*deltaz_ss(i)
        rad_hist_cs = rad_hist_cs + s_hist(i)*rad_ss(i)

        pcld_hc_cs = pcld_hc_cs + s_hc(i)*pcld_ss(i)
        tcld_hc_cs = tcld_hc_cs + s_hc(i)*tcld_ss(i)

        pcld_Cb_cs = pcld_Cb_cs + sCb(i)*pcld_ss(i)
        tcld_Cb_cs = tcld_Cb_cs + sCb(i)*tcld_ss(i)
        em_Cb_cs = em_Cb_cs + sCb(i)*emis_ss(i)

        pcld_ThCi_cs = pcld_ThCi_cs + sThCi(i)*pcld_ss(i)
        tcld_ThCi_cs = tcld_ThCi_cs + sThCi(i)*tcld_ss(i)
        em_ThCi_cs = em_ThCi_cs + sThCi(i)*emis_ss(i)

        pcld_Anv_cs = pcld_Anv_cs + sAnv(i)*pcld_ss(i)
        tcld_Anv_cs = tcld_Anv_cs + sAnv(i)*tcld_ss(i)
        em_Anv_cs = em_Anv_cs + sAnv(i)*emis_ss(i)

        enddo

        deallocate(s)
        deallocate (s_hc)
        deallocate (s_hist)
        deallocate (rneb_max)
        deallocate (emis_ss)
        deallocate (pcld_ss)
        deallocate (tcld_ss)
        deallocate (iwp_ss)
        deallocate (deltaz_ss)
        deallocate (rad_ss)
        deallocate (sCb)
        deallocate (sThCi)
        deallocate (sAnv)

       call normal_undef(pcld_hc_cs,cc_hc_cs)
       call normal_undef(tcld_hc_cs,cc_hc_cs)
       call normal_undef(iwp_hc_cs,cc_hc_cs)
       call normal_undef(em_hc_cs,cc_hc_cs)
       call normal_undef(deltaz_hc_cs,cc_hc_cs)
       
       call normal_undef(iwp_hist_cs,cc_hist_cs)
       call normal_undef(em_hist_cs,cc_hist_cs)
       call normal_undef(deltaz_hist_cs,cc_hist_cs)
       call normal_undef(rad_hist_cs,cc_hist_cs)

       call normal_undef(pcld_Cb_cs,cc_Cb_cs)
       call normal_undef(tcld_Cb_cs,cc_Cb_cs)
       call normal_undef(em_Cb_cs,cc_Cb_cs)

       call normal_undef(pcld_ThCi_cs,cc_ThCi_cs)
       call normal_undef(tcld_ThCi_cs,cc_ThCi_cs)
       call normal_undef(em_ThCi_cs,cc_ThCi_cs)
    
       call normal_undef(pcld_Anv_cs,cc_Anv_cs)
       call normal_undef(tcld_Anv_cs,cc_Anv_cs)
       call normal_undef(em_Anv_cs,cc_Anv_cs)


! Tests

        if (cc_tot_cs > maxval(rneb_cs) .and. &
   abs(cc_tot_cs-maxval(rneb_cs)) > 1.e-4 )  then
          WRITE(abort_message,*) 'cc_tot_cs > max rneb_cs', cc_tot_cs, maxval(rneb_cs)
          CALL abort_physic(modname,abort_message,1)
        endif

        if (iwp_hc_cs < 0.) then
          abort_message= 'cloud_structure: iwp_hc_cs < 0'
          CALL abort_physic(modname,abort_message,1)
        endif
 
        end subroutine cloud_structure


        subroutine normal_undef(num, den)

        REAL, intent(in) :: den
        REAL, intent(inout) :: num

        if (den /= 0) then
        num = num/den
        else
        num = undef
        endif

        end subroutine normal_undef


        subroutine normal2_undef(res,num,den)

        REAL, intent(in) :: den
        REAL, intent(in) :: num
        REAL, intent(out) :: res

        if (den /= 0.) then
        res = num/den
        else
        res = undef
        endif

        end subroutine normal2_undef


        subroutine sous_section(len_cs, rneb_cs, temp_cs, &
   emis_cs, iwco_cs, &
   pres_cs, dz_cs, rhodz_cs, rad_cs, rneb_ord, &
   rnebmax, stot, shc, shist, &
   sCb, sThCi, sAnv, &
   emis, pcld, tcld, iwp, deltaz, rad)

        INTEGER, intent(in) :: len_cs
        REAL, DIMENSION(len_cs), intent(in) :: rneb_cs, temp_cs
        REAL, DIMENSION(len_cs), intent(in) :: emis_cs, iwco_cs, &
   rneb_ord
        REAL, DIMENSION(len_cs), intent(in) :: pres_cs, dz_cs, rad_cs
        REAL, DIMENSION(len_cs), intent(in) :: rhodz_cs
        REAL, DIMENSION(len_cs) :: tau_cs, w
        REAL, intent(in) :: rnebmax
        REAL, intent(inout) :: stot, shc, shist
        REAL, intent(inout) :: sCb, sThCi, sAnv
        REAL, intent(out) :: emis, pcld, tcld, iwp, deltaz, rad

        INTEGER :: i, ideb, ibeg, iend, nuage, visible
        REAL :: som, som_tau, som_iwc, som_dz, som_rad, tau

        CHARACTER (len = 50)      :: modname = 'simu_airs.sous_section'
        CHARACTER (len = 160)     :: abort_message


! Ponderation: 1 pour les nuages, 0 pour les trous

        do i = 1, len_cs
        if (rneb_cs(i) >= rnebmax) then
        w(i) = 1.
        else
        w(i) = 0.
        endif
        enddo

! Calcul des epaisseurs optiques a partir des emissivites

        som = 0.
        do i = 1, len_cs
        if (emis_cs(i) == 1.) then
        tau_cs(i) = 10.
        else
        tau_cs(i) = -log(1.-emis_cs(i))
        endif
        som = som+tau_cs(i)
        enddo


        ideb = 1
        nuage = 0
        visible = 0


! Boucle sur les nuages
        do while (ideb /= 0 .and. ideb <= len_cs)


! Definition d'un nuage de la sous-section

        call topbot(ideb, w, ibeg, iend)
        ideb = iend+1

        if (ibeg > 0) then

        nuage = nuage + 1

! On determine les caracteristiques du nuage
! (ep. optique, ice water path, pression, temperature)

        call caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
   pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
   som_tau, som_iwc, som_dz, som_rad)

! On masque le nuage s'il n'est pas detectable

        call masque(ibeg, iend, som_tau, visible, w)

        endif

! Fin boucle nuages
        enddo


! Analyse du nuage detecte

        call topbot(1, w, ibeg, iend)

        if (ibeg > 0) then

        call caract(ibeg, iend, temp_cs, tau_cs, iwco_cs, &
   pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
   som_tau, som_iwc, som_dz, som_rad)

        tau = som_tau
        emis = 1. - exp(-tau)
        iwp = som_iwc
        deltaz = som_dz
        rad = som_rad

        if (pcld > p_thresh) then

        shc = 0.
        shist = 0.
        sCb = 0.
        sThCi = 0.
        sAnv = 0.

        else

        if (emis < em_min .or. emis > em_max  &
   .or. tcld > 230.) then
        shist = 0.
        endif

        if (emis < 0.98) then
        sCb = 0.
        endif

        if (emis > 0.5 .or. emis < 0.1) then
        sThCi = 0.
        endif

        if (emis <= 0.5 .or. emis >= 0.98) then
        sAnv = 0.
        endif

        endif

        else

        tau = 0.
        emis = 0.
        iwp = 0.
        deltaz = 0.
        pcld = 0.
        tcld = 0.
        stot = 0.
        shc = 0.
        shist = 0.
        rad = 0.
        sCb = 0.
        sThCi = 0.
        sAnv = 0.

        endif


! Tests

        if (iwp < 0.) then
          WRITE(abort_message,*) 'ideb iwp =', ideb, iwp
          CALL abort_physic(modname,abort_message,1)
        endif

        if (deltaz < 0.) then
          WRITE(abort_message,*)'ideb deltaz =', ideb, deltaz
          CALL abort_physic(modname,abort_message,1)
        endif

        if (emis < 0.048 .and. emis /= 0.) then
          WRITE(abort_message,*) 'ideb emis =', ideb, emis
          CALL abort_physic(modname,abort_message,1)
        endif

        end subroutine sous_section


        subroutine masque (ibeg, iend, som_tau, &
   visible, w)

        INTEGER, intent(in) :: ibeg, iend
        REAL, intent(in) :: som_tau

        INTEGER, intent(inout) :: visible
        REAL, DIMENSION(:), intent(inout) :: w

        INTEGER :: i



! Masque

! Cas ou il n'y a pas de nuage visible au-dessus

        if (visible == 0) then

        if (som_tau < tau_thresh) then
        do i = ibeg, iend
        w(i) = 0.
        enddo
        else
        visible = 1
        endif

! Cas ou il y a un nuage visible au-dessus

        else

        do i = ibeg, iend
        w(i) = 0.
        enddo

        endif


        end subroutine masque


         subroutine caract (ibeg, iend, temp_cs, tau_cs, iwco_cs, &
   pres_cs, dz_cs, rhodz_cs, rad_cs, pcld, tcld, &
   som_tau, som_iwc, som_dz, som_rad)

        INTEGER, intent(in) :: ibeg, iend
        REAL, DIMENSION(:), intent(in) :: tau_cs, iwco_cs, temp_cs
        REAL, DIMENSION(:), intent(in) :: pres_cs, dz_cs, rad_cs
        REAL, DIMENSION(:), intent(in) :: rhodz_cs
        REAL, intent(out) :: som_tau, som_iwc, som_dz, som_rad
        REAL , intent(out) :: pcld, tcld

        INTEGER :: i, ibase, imid

        CHARACTER (len = 50)      :: modname = 'simu_airs.caract'
        CHARACTER (len = 160)     :: abort_message

! Somme des epaisseurs optiques et des contenus en glace sur le nuage

        som_tau = 0.
        som_iwc = 0.
        som_dz = 0.
        som_rad = 0.
        ibase = -100

        do i = ibeg, iend

        som_tau = som_tau + tau_cs(i)

        som_dz = som_dz + dz_cs(i)
        som_iwc = som_iwc + iwco_cs(i)*1000*rhodz_cs(i)  ! en g/m2
        som_rad = som_rad + rad_cs(i)*dz_cs(i)

        if (som_tau > 3. .and. ibase == -100) then
        ibase = i-1
        endif

        enddo

        if (som_dz /= 0.) then
          som_rad = som_rad/som_dz
        else
          write(*,*) 'som_dez = 0 STOP'
          write(*,*) 'ibeg, iend =', ibeg, iend
          do i = ibeg, iend
             write(*,*) dz_cs(i), rhodz_cs(i)
          enddo
          abort_message='see above'
          CALL abort_physic(modname,abort_message,1)
        endif

! Determination de Pcld et Tcld

       if (ibase < ibeg) then
       ibase = ibeg
       endif

       imid = (ibeg+ibase)/2

       pcld = pres_cs(imid)/100.        ! pcld en hPa
       tcld = temp_cs(imid)


       end subroutine caract
 
        subroutine topbot(ideb,w,ibeg,iend)

        INTEGER, intent(in) :: ideb
        REAL, DIMENSION(:), intent(in) :: w
        INTEGER, intent(out) :: ibeg, iend

        INTEGER :: i, itest

        itest = 0
        ibeg = 0
        iend = 0

        do i = ideb, size(w)

        if (w(i) == 1. .and. itest == 0) then
        ibeg = i
        itest = 1
        endif

        enddo


        i = ibeg
        do while (w(i) == 1. .and. i <= size(w))
        i = i+1
        enddo
        iend = i-1


        end subroutine topbot

        subroutine ordonne(len_cs, rneb_cs, rneb_ord)

        INTEGER, intent(in) :: len_cs
        REAL, DIMENSION(:), intent(in) :: rneb_cs
        REAL, DIMENSION(:), intent(out) :: rneb_ord

        INTEGER :: i, j, ind_min

        REAL, DIMENSION(len_cs) :: rneb
        REAL :: rneb_max


        do i = 1, size(rneb_cs)
        rneb(i) = rneb_cs(i)
        enddo

        do j = 1, size(rneb_cs)

        rneb_max = 100.

        do i = 1, size(rneb_cs)
        if (rneb(i) < rneb_max) then
        rneb_max = rneb(i)
        ind_min = i
        endif
        enddo

        rneb_ord(j) = rneb_max
        rneb(ind_min) = 100.

        enddo
        
        end subroutine ordonne

 
        subroutine sim_mesh(rneb_1D, temp_1D, emis_1D, &
   iwcon_1D, rad_1D, &
   pres, dz, &
   rhodz_1D, cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, pcld_hc_mesh,&
   tcld_hc_mesh, &
   em_hc_mesh, iwp_hc_mesh, deltaz_hc_mesh, &
   cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
   pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
   pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
   pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
   em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)

       USE dimphy

       REAL, DIMENSION(klev), intent(in) :: rneb_1D, temp_1D, emis_1D, &
   iwcon_1D, rad_1D
        REAL, DIMENSION(klev), intent(in) :: pres, dz, rhodz_1D
        REAL, intent(out) :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
        REAL, intent(out) :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
        REAL, intent(out) :: em_hc_mesh, pcld_hc_mesh, tcld_hc_mesh, &
   iwp_hc_mesh

        REAL, intent(out) :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
        REAL, intent(out) :: pcld_ThCi_mesh, tcld_ThCi_mesh, &
   em_ThCi_mesh
        REAL, intent(out) :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh

        REAL, intent(out) :: em_hist_mesh, iwp_hist_mesh, rad_hist_mesh
        REAL, intent(out) :: deltaz_hc_mesh, deltaz_hist_mesh

        REAL, DIMENSION(:), allocatable :: rneb_cs, temp_cs, emis_cs, &
   iwco_cs
        REAL, DIMENSION(:), allocatable :: pres_cs, dz_cs, rad_cs, &
   rhodz_cs

        INTEGER :: i,j,l
        INTEGER :: ltop, itop, ibot, num_cs, N_cs, len_cs, ics

        REAL :: som_emi_hc,som_pcl_hc,som_tcl_hc,som_iwp_hc,som_hc,&
   som_hist
        REAL :: som_emi_hist, som_iwp_hist, som_deltaz_hc, &
   som_deltaz_hist
        REAL :: som_rad_hist
        REAL :: som_Cb, som_ThCi, som_Anv
        REAL :: som_emi_Cb, som_tcld_Cb, som_pcld_Cb
        REAL :: som_emi_Anv, som_tcld_Anv, som_pcld_Anv
        REAL :: som_emi_ThCi, som_tcld_ThCi, som_pcld_ThCi
        REAL :: tsom_tot, tsom_hc, tsom_hist
        REAL :: prod, prod_hh

        REAL :: cc_tot_cs, cc_hc_cs, cc_hist_cs
        REAL :: cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs
        REAL :: pcld_hc_cs, tcld_hc_cs
        REAL :: em_hc_cs, iwp_hc_cs, deltaz_hc_cs
        REAL :: pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs
        REAL :: pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs
        REAL :: pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs
        REAL :: em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs

        REAL, DIMENSION(klev) :: test_tot, test_hc, test_hist
        REAL, DIMENSION(klev) :: test_pcld, test_tcld, test_em, test_iwp

        CHARACTER (len = 50)      :: modname = 'simu_airs.sim_mesh'
        CHARACTER (len = 160)     :: abort_message
        

        do j = 1, klev
          WRITE(lunout,*) 'simu_airs, j, rneb_1D =', rneb_1D(j)
        enddo

! Definition des structures nuageuses, de la plus haute a la plus basse

        num_cs = 0
        ltop = klev-1

        prod = 1.

        som_emi_hc = 0.
        som_emi_hist = 0.
        som_pcl_hc = 0.
        som_tcl_hc = 0.
        som_iwp_hc = 0.
        som_iwp_hist = 0.
        som_deltaz_hc = 0.
        som_deltaz_hist = 0.
        som_rad_hist = 0.
        som_hc = 0.
        som_hist = 0.

        som_Cb = 0.
        som_ThCi = 0.
        som_Anv = 0.

        som_emi_Cb = 0.
        som_pcld_Cb = 0.
        som_tcld_Cb = 0.

        som_emi_ThCi = 0.
        som_pcld_ThCi = 0.
        som_tcld_ThCi = 0.

        som_emi_Anv = 0.
        som_pcld_Anv = 0.
        som_tcld_Anv = 0.

        tsom_tot = 0.
        tsom_hc = 0.
        tsom_hist = 0.

        do while (ltop >= 1)   ! Boucle sur toute la colonne

        itop = 0

        do l = ltop,1,-1

        if (itop == 0 .and. rneb_1D(l) > 0.001 ) then
        itop = l
        endif

        enddo

        ibot = itop

        do while (rneb_1D(ibot) > 0.001 .and. ibot >= 1)
        ibot = ibot -1
        enddo


        ibot = ibot+1

        if (itop > 0) then    ! itop > 0

        num_cs = num_cs +1
        len_cs = itop-ibot+1

! Allocation et definition des variables de la structure nuageuse
! le premier indice denote ce qui est le plus haut

        allocate (rneb_cs(len_cs))
        allocate (temp_cs(len_cs))
        allocate (emis_cs(len_cs))
        allocate (iwco_cs(len_cs))
        allocate (pres_cs(len_cs))
        allocate (dz_cs(len_cs))
        allocate (rad_cs(len_cs))
        allocate (rhodz_cs(len_cs))

        ics = 0

        do i = itop, ibot, -1
        ics = ics + 1
        rneb_cs(ics) = rneb_1D(i)
        temp_cs(ics) = temp_1D(i)
        emis_cs(ics) = emis_1D(i)
        iwco_cs(ics) = iwcon_1D(i)
        rad_cs(ics) = rad_1D(i)
        pres_cs(ics) = pres(i)
        dz_cs(ics) = dz(i)
        rhodz_cs(ics) = rhodz_1D(i)
        enddo

! Appel du sous_programme cloud_structure

        call cloud_structure(len_cs,rneb_cs,temp_cs,emis_cs,iwco_cs,&
   pres_cs, dz_cs, rhodz_cs, rad_cs, &
   cc_tot_cs, cc_hc_cs, cc_hist_cs, &
   cc_Cb_cs, cc_ThCi_cs, cc_Anv_cs, &
   pcld_hc_cs, tcld_hc_cs, &
   em_hc_cs, iwp_hc_cs, deltaz_hc_cs, &
   pcld_Cb_cs, tcld_Cb_cs, em_Cb_cs, &
   pcld_ThCi_cs, tcld_ThCi_cs, em_ThCi_cs, &
   pcld_Anv_cs, tcld_Anv_cs, em_Anv_cs, &
   em_hist_cs, iwp_hist_cs, deltaz_hist_cs, rad_hist_cs)


        deallocate (rneb_cs)
        deallocate (temp_cs)
        deallocate (emis_cs)
        deallocate (iwco_cs)
        deallocate (pres_cs)
        deallocate (dz_cs)
        deallocate (rad_cs)
        deallocate (rhodz_cs)


! Pour la couverture nuageuse sur la maille

        prod_hh = prod

        prod = prod*(1.-cc_tot_cs)

! Pour les autres variables definies sur la maille

        som_emi_hc = som_emi_hc + em_hc_cs*cc_hc_cs*prod_hh
        som_iwp_hc = som_iwp_hc + iwp_hc_cs*cc_hc_cs*prod_hh
        som_deltaz_hc = som_deltaz_hc + deltaz_hc_cs*cc_hc_cs*prod_hh

        som_emi_Cb = som_emi_Cb + em_Cb_cs*cc_Cb_cs*prod_hh
        som_tcld_Cb = som_tcld_Cb + tcld_Cb_cs*cc_Cb_cs*prod_hh
        som_pcld_Cb = som_pcld_Cb + pcld_Cb_cs*cc_Cb_cs*prod_hh

        som_emi_ThCi = som_emi_ThCi + em_ThCi_cs*cc_ThCi_cs*prod_hh
        som_tcld_ThCi = som_tcld_ThCi + tcld_ThCi_cs*cc_ThCi_cs*prod_hh
        som_pcld_ThCi = som_pcld_ThCi + pcld_ThCi_cs*cc_ThCi_cs*prod_hh

        som_emi_Anv = som_emi_Anv + em_Anv_cs*cc_Anv_cs*prod_hh
        som_tcld_Anv = som_tcld_Anv + tcld_Anv_cs*cc_Anv_cs*prod_hh
        som_pcld_Anv = som_pcld_Anv + pcld_Anv_cs*cc_Anv_cs*prod_hh

        som_emi_hist = som_emi_hist + em_hist_cs*cc_hist_cs*prod_hh
        som_iwp_hist = som_iwp_hist + iwp_hist_cs*cc_hist_cs*prod_hh
        som_deltaz_hist = som_deltaz_hist + &
   deltaz_hist_cs*cc_hist_cs*prod_hh
        som_rad_hist = som_rad_hist + rad_hist_cs*cc_hist_cs*prod_hh

        som_pcl_hc = som_pcl_hc + pcld_hc_cs*cc_hc_cs*prod_hh
        som_tcl_hc = som_tcl_hc + tcld_hc_cs*cc_hc_cs*prod_hh

        som_hc = som_hc + cc_hc_cs*prod_hh
        som_hist = som_hist + cc_hist_cs*prod_hh

        som_Cb = som_Cb + cc_Cb_cs*prod_hh
        som_ThCi = som_ThCi + cc_ThCi_cs*prod_hh
        som_Anv = som_Anv + cc_Anv_cs*prod_hh


! Pour test
        
        call test_bornes('cc_tot_cs     ',cc_tot_cs,1.,0.)
        call test_bornes('cc_hc_cs      ',cc_hc_cs,1.,0.)
        call test_bornes('cc_hist_cs    ',cc_hist_cs,1.,0.)
        call test_bornes('pcld_hc_cs    ',pcld_hc_cs,1200.,0.)
        call test_bornes('tcld_hc_cs    ',tcld_hc_cs,1000.,100.)
        call test_bornes('em_hc_cs      ',em_hc_cs,1000.,0.048)

        test_tot(num_cs) = cc_tot_cs
        test_hc(num_cs) = cc_hc_cs
        test_hist(num_cs) = cc_hist_cs
        test_pcld(num_cs) = pcld_hc_cs
        test_tcld(num_cs) = tcld_hc_cs
        test_em(num_cs) = em_hc_cs
        test_iwp(num_cs) = iwp_hc_cs

        tsom_tot = tsom_tot + cc_tot_cs
        tsom_hc = tsom_hc + cc_hc_cs
        tsom_hist = tsom_hist + cc_hist_cs


        endif                   ! itop > 0

        ltop = ibot -1

        enddo                   ! fin de la boucle sur la colonne

        N_CS = num_cs


! Determination des variables de sortie

        if (N_CS > 0) then   ! if N_CS>0

        cc_tot_mesh = 1. - prod

        cc_hc_mesh = som_hc
        cc_hist_mesh = som_hist

        cc_Cb_mesh = som_Cb
        cc_ThCi_mesh = som_ThCi
        cc_Anv_mesh = som_Anv

        call normal2_undef(pcld_hc_mesh,som_pcl_hc, &
   cc_hc_mesh)
        call normal2_undef(tcld_hc_mesh,som_tcl_hc, &
   cc_hc_mesh)
        call normal2_undef(em_hc_mesh,som_emi_hc, &
   cc_hc_mesh)
        call normal2_undef(iwp_hc_mesh,som_iwp_hc, &
   cc_hc_mesh)
        call normal2_undef(deltaz_hc_mesh,som_deltaz_hc, &
   cc_hc_mesh)

        call normal2_undef(em_Cb_mesh,som_emi_Cb, &
   cc_Cb_mesh)
        call normal2_undef(tcld_Cb_mesh,som_tcld_Cb, &
   cc_Cb_mesh)
        call normal2_undef(pcld_Cb_mesh,som_pcld_Cb, &
   cc_Cb_mesh)

        call normal2_undef(em_ThCi_mesh,som_emi_ThCi, &
   cc_ThCi_mesh)
        call normal2_undef(tcld_ThCi_mesh,som_tcld_ThCi, &
   cc_ThCi_mesh)
        call normal2_undef(pcld_ThCi_mesh,som_pcld_ThCi, &
   cc_ThCi_mesh)

       call normal2_undef(em_Anv_mesh,som_emi_Anv, &
   cc_Anv_mesh)
        call normal2_undef(tcld_Anv_mesh,som_tcld_Anv, &
   cc_Anv_mesh)
        call normal2_undef(pcld_Anv_mesh,som_pcld_Anv, &
   cc_Anv_mesh)


        call normal2_undef(em_hist_mesh,som_emi_hist, &
   cc_hist_mesh)
        call normal2_undef(iwp_hist_mesh,som_iwp_hist, &
   cc_hist_mesh)
        call normal2_undef(deltaz_hist_mesh,som_deltaz_hist, &
   cc_hist_mesh)
        call normal2_undef(rad_hist_mesh,som_rad_hist, &
   cc_hist_mesh)


! Tests 

        ! Tests

       if (cc_tot_mesh > tsom_tot .and. &
   abs(cc_tot_mesh-tsom_tot) > 1.e-4) then
           WRITE(abort_message,*)'cc_tot_mesh > tsom_tot', cc_tot_mesh, tsom_tot
           CALL abort_physic(modname,abort_message,1)
        endif

        if (cc_tot_mesh < maxval(test_tot(1:N_CS)) .and. &
   abs(cc_tot_mesh-maxval(test_tot(1:N_CS))) > 1.e-4) then
           WRITE(abort_message,*) 'cc_tot_mesh < max', cc_tot_mesh, maxval(test_tot(1:N_CS))
           CALL abort_physic(modname,abort_message,1)
        endif

        if (cc_hc_mesh > tsom_hc .and. &
   abs(cc_hc_mesh-tsom_hc) > 1.e-4) then
           WRITE(abort_message,*) 'cc_hc_mesh > tsom_hc', cc_hc_mesh, tsom_hc
           CALL abort_physic(modname,abort_message,1)
        endif

        if (cc_hc_mesh < maxval(test_hc(1:N_CS)) .and. &
   abs(cc_hc_mesh-maxval(test_hc(1:N_CS))) > 1.e-4) then
           WRITE(abort_message,*) 'cc_hc_mesh < max', cc_hc_mesh, maxval(test_hc(1:N_CS))
           CALL abort_physic(modname,abort_message,1)
        endif

        if (cc_hist_mesh > tsom_hist .and. &
   abs(cc_hist_mesh-tsom_hist) > 1.e-4) then
           WRITE(abort_message,*) 'cc_hist_mesh > tsom_hist', cc_hist_mesh, tsom_hist
           CALL abort_physic(modname,abort_message,1)
        endif

        if (cc_hist_mesh < 0.) then
           WRITE(abort_message,*) 'cc_hist_mesh < 0', cc_hist_mesh
           CALL abort_physic(modname,abort_message,1)
        endif

        if ((pcld_hc_mesh > maxval(test_pcld(1:N_CS)) .or. &
   pcld_hc_mesh < minval(test_pcld(1:N_CS))) .and. &
   abs(pcld_hc_mesh-maxval(test_pcld(1:N_CS))) > 1. .and. &
   maxval(test_pcld(1:N_CS)) /= 999. &
   .and. minval(test_pcld(1:N_CS)) /= 999.) then
           WRITE(abort_message,*) 'pcld_hc_mesh est faux', pcld_hc_mesh, maxval(test_pcld(1:N_CS)), &
   minval(test_pcld(1:N_CS))
           CALL abort_physic(modname,abort_message,1)
        endif

       if ((tcld_hc_mesh > maxval(test_tcld(1:N_CS)) .or. &
   tcld_hc_mesh < minval(test_tcld(1:N_CS))) .and. &
   abs(tcld_hc_mesh-maxval(test_tcld(1:N_CS))) > 0.1 .and. &
   maxval(test_tcld(1:N_CS)) /= 999. &
   .and. minval(test_tcld(1:N_CS)) /= 999.) then
           WRITE(abort_message,*) 'tcld_hc_mesh est faux', tcld_hc_mesh, maxval(test_tcld(1:N_CS)), &
   minval(test_tcld(1:N_CS))
           CALL abort_physic(modname,abort_message,1)
        endif

        if ((em_hc_mesh > maxval(test_em(1:N_CS)) .or. &
   em_hc_mesh < minval(test_em(1:N_CS))) .and. &
   abs(em_hc_mesh-maxval(test_em(1:N_CS))) > 1.e-4 .and. &
   minval(test_em(1:N_CS)) /= 999. .and. &
   maxval(test_em(1:N_CS)) /= 999. ) then
           WRITE(abort_message,*) 'em_hc_mesh est faux', em_hc_mesh, maxval(test_em(1:N_CS)), &
   minval(test_em(1:N_CS))
           CALL abort_physic(modname,abort_message,1)
        endif

        else               ! if N_CS>0

        cc_tot_mesh = 0.
        cc_hc_mesh = 0.
        cc_hist_mesh = 0.

        cc_Cb_mesh = 0.
        cc_ThCi_mesh = 0.
        cc_Anv_mesh = 0.

        iwp_hc_mesh = undef 
        deltaz_hc_mesh = undef 
        em_hc_mesh = undef 
        iwp_hist_mesh = undef 
        deltaz_hist_mesh = undef 
        rad_hist_mesh = undef 
        em_hist_mesh = undef 
        pcld_hc_mesh = undef 
        tcld_hc_mesh = undef 

        pcld_Cb_mesh = undef 
        tcld_Cb_mesh = undef 
        em_Cb_mesh = undef 

        pcld_ThCi_mesh = undef 
        tcld_ThCi_mesh = undef 
        em_ThCi_mesh = undef 

        pcld_Anv_mesh = undef 
        tcld_Anv_mesh = undef 
        em_Anv_mesh = undef 


        endif                  ! if N_CS>0

        end subroutine sim_mesh

        subroutine test_bornes(sx,x,bsup,binf)

        REAL, intent(in) :: x, bsup, binf
        character*14, intent(in) :: sx
        CHARACTER (len = 50)      :: modname = 'simu_airs.test_bornes'
        CHARACTER (len = 160)     :: abort_message

        if (x > bsup .or. x < binf) then
          WRITE(abort_message,*) sx, 'est faux', sx, x
          CALL abort_physic(modname,abort_message,1)
        endif
 
        end subroutine test_bornes

        end module m_simu_airs


        subroutine simu_airs &
   (itap, rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, rad_airs, &
   geop_airs, pplay_airs, paprs_airs, &
   map_prop_hc,map_prop_hist,&
   map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
   map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
   map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
   map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
   map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
   map_ntot,map_hc,map_hist,&
   map_Cb,map_ThCi,map_Anv,alt_tropo )

        USE dimphy
        USE m_simu_airs

        IMPLICIT NONE

        include "YOMCST.h"

        INTEGER,intent(in) :: itap

        REAL, DIMENSION(klon,klev), intent(in) :: &
   rneb_airs, temp_airs, cldemi_airs, iwcon0_airs, &
   rad_airs, geop_airs, pplay_airs, paprs_airs

       REAL, DIMENSION(klon,klev) :: &
   rhodz_airs, rho_airs, iwcon_airs

        REAL, DIMENSION(klon),intent(out) :: alt_tropo

        REAL, DIMENSION(klev) :: rneb_1D, temp_1D, &
   emis_1D, rad_1D, pres_1D, alt_1D, &
   rhodz_1D, dz_1D, iwcon_1D

        INTEGER :: i, j

        REAL :: cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
        REAL :: cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh
        REAL :: pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
        REAL :: em_hist_mesh, iwp_hist_mesh
        REAL :: deltaz_hc_mesh, deltaz_hist_mesh, rad_hist_mesh
        REAL :: pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh
        REAL :: pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh
        REAL :: pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh

        REAL, DIMENSION(klon),intent(out) :: map_prop_hc, map_prop_hist
        REAL, DIMENSION(klon),intent(out) :: map_emis_hc, map_iwp_hc
        REAL, DIMENSION(klon),intent(out) :: map_deltaz_hc, map_pcld_hc
        REAL, DIMENSION(klon),intent(out) :: map_tcld_hc
        REAL, DIMENSION(klon),intent(out) :: map_emis_Cb,map_pcld_Cb,map_tcld_Cb 
        REAL, DIMENSION(klon),intent(out) :: &
   map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi
        REAL, DIMENSION(klon),intent(out) :: &
   map_emis_Anv,map_pcld_Anv,map_tcld_Anv
        REAL, DIMENSION(klon),intent(out) :: &
   map_emis_hist,map_iwp_hist,map_deltaz_hist,&
   map_rad_hist
        REAL, DIMENSION(klon),intent(out) :: map_ntot,map_hc,map_hist
        REAL, DIMENSION(klon),intent(out) :: map_Cb,map_ThCi,map_Anv
 
 
        write(*,*) 'simu_airs'
        write(*,*) 'itap, klon, klev', itap, klon, klev
        write(*,*) 'RG, RD =', RG, RD


! Definition des variables 1D

        do i = 1, klon
        do j = 1, klev-1
        rhodz_airs(i,j) = &
   (paprs_airs(i,j)-paprs_airs(i,j+1))/RG
        enddo
        rhodz_airs(i,klev) = 0.
        enddo

        do i = 1, klon
        do j = 1,klev
        rho_airs(i,j) = &
   pplay_airs(i,j)/(temp_airs(i,j)*RD)

        if (rneb_airs(i,j) > 0.001) then
        iwcon_airs(i,j) = iwcon0_airs(i,j)/rneb_airs(i,j)
        else
        iwcon_airs(i,j) = 0.
        endif
 
        enddo
        enddo

!=============================================================================

        do i = 1, klon  ! boucle sur les points de grille

!=============================================================================
        
        do j = 1,klev

        rneb_1D(j) = rneb_airs(i,j)
        temp_1D(j) = temp_airs(i,j)
        emis_1D(j) = cldemi_airs(i,j)
        iwcon_1D(j) = iwcon_airs(i,j)
        rad_1D(j) = rad_airs(i,j)
        pres_1D(j) = pplay_airs(i,j)
        alt_1D(j) = geop_airs(i,j)/RG
        rhodz_1D(j) = rhodz_airs(i,j)
        dz_1D(j) = rhodz_airs(i,j)/rho_airs(i,j)

        enddo

        alt_tropo(i) = &
   search_tropopause(pres_1D/100.,temp_1D,alt_1D,klev)


! Appel du ss-programme sim_mesh

!        if (itap .eq. 1 ) then

        call sim_mesh(rneb_1D, temp_1D, emis_1D, iwcon_1D, rad_1D, &
   pres_1D, dz_1D, rhodz_1D, &
   cc_tot_mesh, cc_hc_mesh, cc_hist_mesh, &
   pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh, &
   deltaz_hc_mesh,&
   cc_Cb_mesh, cc_ThCi_mesh, cc_Anv_mesh, &
   pcld_Cb_mesh, tcld_Cb_mesh, em_Cb_mesh, &
   pcld_ThCi_mesh, tcld_ThCi_mesh, em_ThCi_mesh, &
   pcld_Anv_mesh, tcld_Anv_mesh, em_Anv_mesh, &
   em_hist_mesh, iwp_hist_mesh, deltaz_hist_mesh, rad_hist_mesh)

         write(*,*) '===================================='
         write(*,*) 'itap, i:', itap, i 
         write(*,*) 'cc_tot, cc_hc, cc_hist, pcld_hc, tcld_hc, em_hc, &
   iwp_hc, em_hist, iwp_hist ='
         write(*,*) cc_tot_mesh, cc_hc_mesh, cc_hist_mesh
         write(*,*) pcld_hc_mesh, tcld_hc_mesh, em_hc_mesh, iwp_hc_mesh
         write(*,*)  em_hist_mesh, iwp_hist_mesh

!        endif

! Definition des variables a ecrire dans le fichier de sortie

        call normal2_undef(map_prop_hc(i),cc_hc_mesh, &
   cc_tot_mesh)
        call normal2_undef(map_prop_hist(i),cc_hist_mesh, &
   cc_tot_mesh)

       map_emis_hc(i) = em_hc_mesh
       map_iwp_hc(i) = iwp_hc_mesh
       map_deltaz_hc(i) = deltaz_hc_mesh
       map_pcld_hc(i) = pcld_hc_mesh
       map_tcld_hc(i) = tcld_hc_mesh

       map_emis_Cb(i) = em_Cb_mesh
       map_pcld_Cb(i) = pcld_Cb_mesh
       map_tcld_Cb(i) = tcld_Cb_mesh

       map_emis_ThCi(i) = em_ThCi_mesh
       map_pcld_ThCi(i) = pcld_ThCi_mesh
       map_tcld_ThCi(i) = tcld_ThCi_mesh

       map_emis_Anv(i) = em_Anv_mesh
       map_pcld_Anv(i) = pcld_Anv_mesh
       map_tcld_Anv(i) = tcld_Anv_mesh

       map_emis_hist(i) = em_hist_mesh
       map_iwp_hist(i) = iwp_hist_mesh
       map_deltaz_hist(i) = deltaz_hist_mesh
       map_rad_hist(i) = rad_hist_mesh

       map_ntot(i) = cc_tot_mesh
       map_hc(i) = cc_hc_mesh
       map_hist(i) = cc_hist_mesh

       map_Cb(i) = cc_Cb_mesh
       map_ThCi(i) = cc_ThCi_mesh
       map_Anv(i) = cc_Anv_mesh


        enddo         ! fin boucle sur les points de grille

        

        end subroutine simu_airs