c
c===========================================================================
c
c Fabrication des profils initiaux pour les simulations du cas idealise pour
c l'etude de la sensibilite de la convection a l'humidite.
c (S.H. Derbyshire et al., 2004: Sensitivity of moist convection to 
c  environmental humidity.  Q. J. Roy. Meteorol. Soc., 130, 3055-3079).
c
c===========================================================================
c
c    ztop ----------------------------------------------------------------
c              ^                          ^                            .
c              |                          |               .           .
c              |                          |                          .
c              |                          |                .        .theta
c              |                          |           RH (q=cste)  .
c              |                          |                 .     .
c              |                          |                      .
c              |                          |                  .  .
c              |                          |                    .
c              |                          |                   .
c              |                          |                  . 
c    ztop2 -----------------------------------------------------.----------
c              |                 ^        |                 .   .
c              |                 |        |                     .
c          dz =|zstep2           |        |                .    .
c              |                 |        |nblvlm               .
c              |                 |        |               .     .
c              |                 |        |                     .RH=rhzmid
c              |                 |        |       theta  .      .
c              |                 |        |                     .
c              |                 |        |             .       . 
c              |                 |        |                     .
c              |                 |        |            .        .
c              |                 |        |                     .
c              |                 |        |           .         .
c              V                 |        |                     .
c    ztop1  -----------------------------------------.---------------.----
c              ^        ^        |        |                          .
c              |        |        |nblvlm2 |         .                .
c              |        |        |        |       theta              .
c              |        |        |        |        .                 .RH=rhz0
c              |        |        |        |                          .
c    ztopbl --------------------------------------.------------------.-------
c              |        |nblvlm1 |        |       .                  .
c          dz =|zstep1  |        |        |       .                  .
c              |        |        |        |       .                  .
c              |        |        |        |       .                  .
c              V        V        V        V       .                  .
c    z=0  ----------------------------------------------------------------
c         ///////////////////////////////////////////////////////////////
c
      IMPLICIT none
c
c
      REAL val
      REAL plim
      REAL ak
      INTEGER ierr
      INTEGER i
      INTEGER lu
      INTEGER k
      INTEGER mlz
      INTEGER mlzh
      INTEGER j
c
      INTEGER nblvlm1, nblvlm2, nblvlm
      REAL zstep1,zstep2
      REAL tetaz0,tetaztop, tetaztop2
      REAL ztop, ztopbl, ztop1, ztop2
      REAL dztop, dztop2
c
      REAL psolpa
c
      REAL Tsol
      REAL dtetadz2, dtetadz
c
      REAL qsol
      REAL rhz0,rhzmid
c
      REAL rz0,G,CPD,RD,Akappa,fact
      REAL pkappa(100),pkappasol,P0
      REAL playm(100),hplaym(100)
c
      REAL e12
c
      REAL tempz(100),ovapz(100),rhz(100),qlz(100),tetaz(100)
     $           ,uz(100),vz(100),zz(100)

c
c  -----------------------------------------------------
c       Statement function : saturation vapour pressure
c
      REAL qsatl,qsats,ptarg
c       Attention : qsatl,qsats donne 0.622*Psat
      qsats(ptarg) = 100.0 * 0.622 * 10.0
     .           ** (2.07023 - 0.00320991 * ptarg
     .           - 2484.896 / ptarg + 3.56654 * LOG10(ptarg))
      qsatl(ptarg) = 100.0 * 0.622 * 10.0
     .           ** (23.8319 - 2948.964 / ptarg
     .           - 5.028 * LOG10(ptarg)
     .           - 29810.16 * EXP( - 0.0699382 * ptarg)
     .           + 25.21935 * EXP( - 2999.924 / ptarg))
c  -----------------------------------------------------
c
c
       nblvlm1 = 40
       nblvlm2 = 50
       nblvlm  = 70
cc       zstep1 = 60.
cc       zstep2 = 1300.
c
       tetaz0=293.
       dtetadz2 = 3.e-3    ! K/m ; gradient de theta de ztopbl a ztop2
       dtetadz  = 40.e-3   ! K/m ; gradient de theta de ztop2 a ztop
   
       rhz0    = .8          ! humidite relative entre ztopbl et ztop1
       rhzmid  = .80       ! valeurs possibles : .25, .5, .7, .9
       ztopbl  = 1000.     ! sommet couche limite
       ztop1   = 2000.     ! sommet de la couche humide basse de la tropo libre.
       ztop2   = 15000.      ! sommet profil lineaire pour Theta
       ztop    = 40000.     ! sommet grille
c
       open (99,file='profile.input',status='old',iostat=ierr)
        if (ierr /= 0) then
            print *,' Pb ouverture profile.input'
            stop
        else
           read (99,*,iostat=ierr) rhzmid
           rhzmid = rhzmid/100.
           if (ierr /= 0) then
             print *,' Pb lecture profile.input'
             stop
           else
             print *,'Humidite tropo libre = ',rhzmid
           endif 
           close(99)
        endif
c
       zstep1 = ztop1/nblvlm1
       zstep2 = (ztop-ztop1)/(nblvlm-nblvlm1)
       dztop2  = ztop2 - ztopbl
       dztop   = ztop  - ztop2
ccc       tetaztop = 335.
       tetaztop2 = tetaz0    + dtetadz2*dztop2
       tetaztop  = tetaztop2 + dtetadz *dztop
c
       G = 9.81
       RD = 1000.0 * 1.380658E-23 * 6.0221367E+23 / 28.9644
       CPD = 3.5 * RD
       AKAPPA = 2./7
       P0 = 1.e5
       fact = P0**akappa*G*.5/CPD
c
       rz0 = 0.1          ! longueur de rugosite
       psolpa = 1.e5
       Tsol = 294.
       qsol = qsatl(Tsol)/psolpa
c
       pkappasol = P0**akappa
       zz(1) = zstep1/2.
       tetaz(1) = tetaz0
       pkappa(1) = pkappasol - fact*(2./tetaz0)*zz(1)
c
       print *,'zstep1, zstep2 ',zstep1, zstep2
       print *, ' ztop2, ztop ',ztop2, ztop
       print *, 'tetaztop2, tetaztop ',tetaztop2, tetaztop
c
c        Calcul de Theta et de P**kappa
c
      DO k = 2,nblvlm2
        IF (k .le. nblvlm1) THEN
          zz(k) = zz(k-1) + zstep1
        ELSE
          zz(k) = zz(k-1) + zstep2
        ENDIF
        tetaz(k) = tetaz0
     $           + MAX(0.,(zz(k)-ztopbl)*dtetadz2)
        pkappa(k) = pkappa(k-1) 
     $           - fact*(1./tetaz(k-1)+1./tetaz(k))*(zz(k)-zz(k-1))
        print *, 'k,zz(k),tetaz(k) ',k,zz(k),tetaz(k)
      ENDDO
      DO k = nblvlm2+1,nblvlm
        zz(k) = zz(k-1) + zstep2
        tetaz(k) = tetaztop2
     $           + MAX(0.,(zz(k)-ztop2)*dtetadz)
        pkappa(k) = pkappa(k-1) 
     $           - fact*(1./tetaz(k-1)+1./tetaz(k))*(zz(k)-zz(k-1))
        print *, 'k,zz(k),tetaz(k) ',k,zz(k),tetaz(k)
      ENDDO
c
c      Passage aux pressions
c
      DO k = 1,nblvlm2
        IF (zz(k) .lt. ztop1) THEN
          rhz(k) = rhz0
        ELSE
          rhz(k) = rhzmid
        ENDIF
         playm(k) = pkappa(k)**(1./akappa)
         tempz(k) = tetaz(k)*pkappa(k)/P0**akappa
         IF (tempz(k) .lt. 273.15) THEN
           ovapz(k) = rhz(k)*qsats(tempz(k))/playm(k)
         ELSE
           ovapz(k) = rhz(k)*qsatl(tempz(k))/playm(k)
         ENDIF
         uz(k) = 0.5*ALOG(1.+ zz(k)/rz0)
         vz(k) = 0.
      ENDDO
      DO k = nblvlm2+1,nblvlm
        ovapz(k) = ovapz(nblvlm2)
        playm(k) = pkappa(k)**(1./akappa)
        tempz(k) = tetaz(k)*pkappa(k)/P0**akappa
         IF (tempz(k) .lt. 273.15) THEN
           rhz(k) = ovapz(k)*playm(k)/qsats(tempz(k))
         ELSE
           rhz(k) = ovapz(k)*playm(k)/qsatl(tempz(k))
         ENDIF
         uz(k) = 0.5*ALOG(1.+ zz(k)/rz0)
         vz(k) = 0.
       ENDDO

c
       open (99,file='prof.inp',status='unknown',iostat=ierr)
        if (ierr /= 0) then
            print *,' Pb ecriture prof.inp'
            stop
        else
           e12 = 1.
           write (99,*) nblvlm
           do k=1,nblvlm
             write (99,1002) zz(k),tempz(k),ovapz(k),                         &
     &                      uz(k),vz(k),e12
           enddo
           close(99)
        endif

ccc      Si la pression est en HPa, la multiplier par 100
cc      if (playm(1) .lt. 10000.) then
cc        do mlz = 1,nblvlm
cc         playm(mlz) = playm(mlz)*100.
cc        enddo
cc      endif
      print*,(playm(mlz),mlz=1,nblvlm)
c
 1000 format (a4)
 1001 format(5x,i2)
 1002 format (2(2x,f10.4),2x,e13.6,3(2x,f10.5))
c 
      print*,' '
      do mlzh=1,nblvlm
      hplaym(mlzh)=playm(mlzh)/100.
      enddo
c
      print*,'pression en hPa de chaque couche du meso-NH: '
      print*,(hplaym(mlzh),mlzh=1,nblvlm)
c
      return
c
      end

