source: LMDZ6/branches/Amaury_dev/libf/phylmd/dyn1d/lmdz_old_1dconv.f90

MODULE lmdz_old_1dconv
  PRIVATE  ! -- We'd love to put IMPLICIT NONE;  here...
  PUBLIC get_uvd, copie, get_uvd2, rdgrads, spaces

  REAL play(100)  !pression en Pa au milieu de chaque couche GCM
  INTEGER JM(100) !pression en Pa au milieu de chaque couche GCM
  REAL coef1(100) !coefficient d interpolation
  REAL coef2(100) !coefficient d interpolation
  INTEGER klev

  INTEGER nblvlm !nombre de niveau de pression du mesoNH
  REAL playm(100)  !pression en Pa au milieu de chaque couche Meso-NH
  REAL hplaym(100) !pression en hPa milieux des couches Meso-NH


CONTAINS

  subroutine get_uvd(itap, dtime, file_forctl, file_fordat, &
          &       ht, hq, hw, hu, hv, hthturb, hqturb, &
          &       Ts, imp_fcg, ts_fcg, Tp_fcg, Turb_fcg)

    USE lmdz_yomcst

    IMPLICIT NONE

    !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    ! cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
    ! pouvoir calculer la convergence et le cisaillement dans la physiq
    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    INTEGER i, j, k, ll, in
    CHARACTER*80 file_forctl, file_fordat

    !======================================================================
    ! methode: on va chercher les donnees du mesoNH de meteo france, on y
    !          a acces a tout pas detemps grace a la routine rdgrads qui
    !          est une boucle lisant dans ces fichiers.
    !          Puis on interpole ces donnes sur les 11 niveaux du gcm et
    !          et sur les pas de temps de ce meme gcm
    !----------------------------------------------------------------------
    ! input:
    !       pasmax     :nombre de pas de temps maximum du mesoNH
    !       dt         :pas de temps du meso_NH (en secondes)
    !----------------------------------------------------------------------
    INTEGER pasmax, dt
    save pasmax, dt
    !----------------------------------------------------------------------
    ! arguments:
    !           itap   :compteur de la physique(le nombre de ces pas est
    !                   fixe dans la subroutine calcul_ini_gcm de interpo
    !                   -lation
    !           dtime  :pas detemps du gcm (en secondes)
    !           ht     :convergence horizontale de temperature(K/s)
    !           hq     :    "         "       d humidite (kg/kg/s)
    !           hw     :vitesse verticale moyenne (m/s**2)
    !           hu     :convergence horizontale d impulsion le long de x
    !                  (kg/(m^2 s^2)
    !           hv     : idem le long de y.
    !           Ts     : Temperature de surface (K)
    !           imp_fcg: var. LOGICAL .EQ. T si forcage en impulsion
    !           ts_fcg: var. LOGICAL .EQ. T si forcage en Ts present dans fichier
    !           Tp_fcg: var. LOGICAL .EQ. T si forcage donne en Temp potentielle
    !           Turb_fcg: var. LOGICAL .EQ. T si forcage turbulent present dans fichier
    !----------------------------------------------------------------------
    INTEGER itap
    REAL dtime
    REAL ht(:)
    REAL hq(:)
    REAL hu(:)
    REAL hv(:)
    REAL hw(:)
    REAL hthturb(:)
    REAL hqturb(:)
    REAL Ts, Ts_subr
    LOGICAL imp_fcg
    LOGICAL ts_fcg
    LOGICAL Tp_fcg
    LOGICAL Turb_fcg
    !----------------------------------------------------------------------
    ! Variables internes de get_uvd (note : l interpolation temporelle
    ! est faite entre les pas de temps before et after, sur les variables
    ! definies sur la grille du SCM; on atteint exactement les valeurs Meso
    ! aux milieux des pas de temps Meso)
    !     time0     :date initiale en secondes
    !     time      :temps associe a chaque pas du SCM
    !     pas       :numero du pas du meso_NH (on lit en pas : le premier pas
    !                 des donnees est duplique)
    !     pasprev   :numero du pas de lecture precedent
    !     htaft     :advection horizontale de temp. au pas de temps after
    !     hqaft     :    "         "      d humidite        "
    !     hwaft     :vitesse verticalle moyenne  au pas de temps after
    !     huaft,hvaft :advection horizontale d impulsion au pas de temps after
    !     tsaft     : surface temperature 'after time step'
    !     htbef     :idem htaft, mais pour le pas de temps before
    !     hqbef     :voir hqaft
    !     hwbef     :voir hwaft
    !     hubef,hvbef : idem huaft,hvaft, mais pour before
    !     tsbef     : surface temperature 'before time step'
    !----------------------------------------------------------------------
    INTEGER time0, pas, pasprev
    save time0, pas, pasprev
    REAL time
    REAL htaft(100), hqaft(100), hwaft(100), huaft(100), hvaft(100)
    REAL hthturbaft(100), hqturbaft(100)
    REAL Tsaft
    save htaft, hqaft, hwaft, huaft, hvaft, hthturbaft, hqturbaft
    REAL htbef(100), hqbef(100), hwbef(100), hubef(100), hvbef(100)
    REAL hthturbbef(100), hqturbbef(100)
    REAL Tsbef
    save htbef, hqbef, hwbef, hubef, hvbef, hthturbbef, hqturbbef

    REAL timeaft, timebef
    save timeaft, timebef
    INTEGER temps
    CHARACTER*4 string
    !----------------------------------------------------------------------
    ! variables arguments de la subroutine rdgrads
    !---------------------------------------------------------------------
    INTEGER icompt, icomp1 !compteurs de rdgrads
    REAL z(100)         ! altitude (grille Meso)
    REAL ht_mes(100)    !convergence horizontale de temperature
    !-(grille Meso)
    REAL hq_mes(100)    !convergence horizontale d humidite
    !(grille Meso)
    REAL hw_mes(100)    !vitesse verticale moyenne
    !(grille Meso)
    REAL hu_mes(100), hv_mes(100)    !convergence horizontale d impulsion
    !(grille Meso)
    REAL hthturb_mes(100) !tendance horizontale de T_pot, due aux
    !flux turbulents
    REAL hqturb_mes(100) !tendance horizontale d humidite, due aux
    !flux turbulents

    !---------------------------------------------------------------------
    ! variable argument de la subroutine copie
    !---------------------------------------------------------------------
    ! SB        real pplay(100)    !pression en milieu de couche du gcm
    ! SB                            !argument de la physique
    !---------------------------------------------------------------------
    ! variables destinees a la lecture du pas de temps du fichier de donnees
    !---------------------------------------------------------------------
    CHARACTER*80 aaa, atemps, apasmax
    INTEGER nch, imn, ipa
    !---------------------------------------------------------------------
    PRINT*, 'le pas itap est:', itap
    !*** on determine le pas du meso_NH correspondant au nouvel itap ***
    !*** pour aller chercher les champs dans rdgrads                 ***

    time = time0 + itap * dtime
    !c        temps=int(time/dt+1)
    !c        pas=min(temps,pasmax)
    temps = 1 + int((dt + 2 * time) / (2 * dt))
    pas = min(temps, pasmax - 1)
    PRINT*, 'le pas Meso est:', pas


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

    !*** on remplit les champs before avec les champs after du pas   ***
    !*** precedent en format gcm                                     ***
    IF(pas>pasprev)THEN
      DO i = 1, klev
        htbef(i) = htaft(i)
        hqbef(i) = hqaft(i)
        hwbef(i) = hwaft(i)
        hubef(i) = huaft(i)
        hvbef(i) = hvaft(i)
        hThTurbbef(i) = hThTurbaft(i)
        hqTurbbef(i) = hqTurbaft(i)
      enddo
      tsbef = tsaft
      timebef = pasprev * dt
      timeaft = timebef + dt
      icomp1 = nblvlm * 4
      IF (ts_fcg) icomp1 = icomp1 + 1
      IF (imp_fcg) icomp1 = icomp1 + nblvlm * 2
      IF (Turb_fcg) icomp1 = icomp1 + nblvlm * 2
      icompt = icomp1 * pas
      PRINT *, 'imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt'
      PRINT *, imp_fcg, ts_fcg, Turb_fcg, pas, nblvlm, icompt
      PRINT*, 'le pas pas est:', pas
      !*** on va chercher les nouveaux champs after dans toga.dat     ***
      !*** champs en format meso_NH                                   ***
      open(99, FILE = file_fordat, FORM = 'UNFORMATTED', &
              &             ACCESS = 'DIRECT', RECL = 8)
      CALL rdgrads(99, icompt, nblvlm, z, ht_mes, hq_mes, hw_mes                &
              &, hu_mes, hv_mes, hthturb_mes, hqturb_mes                 &
              &, ts_fcg, ts_subr, imp_fcg, Turb_fcg)

      IF(Tp_fcg) THEN
        !     (le forcage est donne en temperature potentielle)
        DO i = 1, nblvlm
          ht_mes(i) = ht_mes(i) * (hplaym(i) / 1000.)**rkappa
        enddo
      endif ! Tp_fcg
      IF(Turb_fcg) THEN
        DO i = 1, nblvlm
          hThTurb_mes(i) = hThTurb_mes(i) * (hplaym(i) / 1000.)**rkappa
        enddo
      endif  ! Turb_fcg

      PRINT*, 'ht_mes ', (ht_mes(i), i = 1, nblvlm)
      PRINT*, 'hq_mes ', (hq_mes(i), i = 1, nblvlm)
      PRINT*, 'hw_mes ', (hw_mes(i), i = 1, nblvlm)
      IF(imp_fcg) THEN
        PRINT*, 'hu_mes ', (hu_mes(i), i = 1, nblvlm)
        PRINT*, 'hv_mes ', (hv_mes(i), i = 1, nblvlm)
      endif
      IF(Turb_fcg) THEN
        PRINT*, 'hThTurb_mes ', (hThTurb_mes(i), i = 1, nblvlm)
        PRINT*, 'hqTurb_mes ', (hqTurb_mes(i), i = 1, nblvlm)
      endif
      IF (ts_fcg) PRINT*, 'ts_subr', ts_subr
      !*** on interpole les champs meso_NH sur les niveaux de pression***
      !*** gcm . on obtient le nouveau champ after                    ***
      DO k = 1, klev
        IF (JM(k) == 0) THEN
          htaft(k) = ht_mes(jm(k) + 1)
          hqaft(k) = hq_mes(jm(k) + 1)
          hwaft(k) = hw_mes(jm(k) + 1)
          IF(imp_fcg) THEN
            huaft(k) = hu_mes(jm(k) + 1)
            hvaft(k) = hv_mes(jm(k) + 1)
          endif ! imp_fcg
          IF(Turb_fcg) THEN
            hThTurbaft(k) = hThTurb_mes(jm(k) + 1)
            hqTurbaft(k) = hqTurb_mes(jm(k) + 1)
          endif ! Turb_fcg
        else ! JM(k) .EQ. 0
          htaft(k) = coef1(k) * ht_mes(jm(k)) + coef2(k) * ht_mes(jm(k) + 1)
          hqaft(k) = coef1(k) * hq_mes(jm(k)) + coef2(k) * hq_mes(jm(k) + 1)
          hwaft(k) = coef1(k) * hw_mes(jm(k)) + coef2(k) * hw_mes(jm(k) + 1)
          IF(imp_fcg) THEN
            huaft(k) = coef1(k) * hu_mes(jm(k)) + coef2(k) * hu_mes(jm(k) + 1)
            hvaft(k) = coef1(k) * hv_mes(jm(k)) + coef2(k) * hv_mes(jm(k) + 1)
          endif ! imp_fcg
          IF(Turb_fcg) THEN
            hThTurbaft(k) = coef1(k) * hThTurb_mes(jm(k))                            &
                    & + coef2(k) * hThTurb_mes(jm(k) + 1)
            hqTurbaft(k) = coef1(k) * hqTurb_mes(jm(k))                             &
                    & + coef2(k) * hqTurb_mes(jm(k) + 1)
          endif ! Turb_fcg
        endif ! JM(k) .EQ. 0
      enddo
      tsaft = ts_subr
      pasprev = pas
    else ! pas.gt.pasprev
      PRINT*, 'timebef est:', timebef
    endif ! pas.gt.pasprev    fin du bloc relatif au passage au pas
    !de temps (meso) suivant
    !*** si on atteint le pas max des donnees experimentales ,on     ***
    !*** on conserve les derniers champs calcules                    ***
    IF(temps>=pasmax)THEN
      DO ll = 1, klev
        ht(ll) = htaft(ll)
        hq(ll) = hqaft(ll)
        hw(ll) = hwaft(ll)
        hu(ll) = huaft(ll)
        hv(ll) = hvaft(ll)
        hThTurb(ll) = hThTurbaft(ll)
        hqTurb(ll) = hqTurbaft(ll)
      enddo
      ts_subr = tsaft
    else ! temps.ge.pasmax
      !*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
      !** des pas de temps de 1h du meso_NH                            ***
      DO j = 1, klev
        ht(j) = ((timeaft - time) * htbef(j) + (time - timebef) * htaft(j)) / dt
        hq(j) = ((timeaft - time) * hqbef(j) + (time - timebef) * hqaft(j)) / dt
        hw(j) = ((timeaft - time) * hwbef(j) + (time - timebef) * hwaft(j)) / dt
        IF(imp_fcg) THEN
          hu(j) = ((timeaft - time) * hubef(j) + (time - timebef) * huaft(j)) / dt
          hv(j) = ((timeaft - time) * hvbef(j) + (time - timebef) * hvaft(j)) / dt
        endif ! imp_fcg
        IF(Turb_fcg) THEN
          hThTurb(j) = ((timeaft - time) * hThTurbbef(j)                           &
                  & + (time - timebef) * hThTurbaft(j)) / dt
          hqTurb(j) = ((timeaft - time) * hqTurbbef(j)                            &
                  & + (time - timebef) * hqTurbaft(j)) / dt
        endif ! Turb_fcg
      enddo
      ts_subr = ((timeaft - time) * tsbef + (time - timebef) * tsaft) / dt
    endif ! temps.ge.pasmax

    PRINT *, ' time,timebef,timeaft', time, timebef, timeaft
    PRINT *, ' ht,htbef,htaft,hthturb,hthturbbef,hthturbaft'
    DO j = 1, klev
      PRINT *, j, ht(j), htbef(j), htaft(j), &
              &             hthturb(j), hthturbbef(j), hthturbaft(j)
    enddo
    PRINT *, ' hq,hqbef,hqaft,hqturb,hqturbbef,hqturbaft'
    DO j = 1, klev
      PRINT *, j, hq(j), hqbef(j), hqaft(j), &
              &             hqturb(j), hqturbbef(j), hqturbaft(j)
    enddo

    !-------------------------------------------------------------------

    IF (Ts_fcg) Ts = Ts_subr
    return

    !-----------------------------------------------------------------------
    ! on sort les champs de "convergence" pour l instant initial 'in'
    ! ceci se passe au pas temps itap=0 de la physique
    !-----------------------------------------------------------------------
    entry get_uvd2(itap, dtime, file_forctl, file_fordat, &
            &           ht, hq, hw, hu, hv, hThTurb, hqTurb, ts, &
            &           imp_fcg, ts_fcg, Tp_fcg, Turb_fcg)
    PRINT*, 'le pas itap est:', itap

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

    WRITE(*, '(a)') 'OPEN ' // file_forctl
    open(97, FILE = file_forctl, FORM = 'FORMATTED')

    !------------------
    DO i = 1, 1000
      read(97, 1000, end = 999) string
      1000 format (a4)
      IF (string == 'TDEF') go to 50
    enddo
    50   backspace(97)
    !-------------------------------------------------------------------
    !   *** on lit le pas de temps dans le fichier de donnees ***
    !   *** "forcing.ctl" et pasmax                           ***
    !-------------------------------------------------------------------
    read(97, 2000) aaa
    2000  format (a80)
    PRINT*, 'aaa est', aaa
    aaa = spaces(aaa, 1)
    PRINT*, 'aaa', aaa
    CALL getsch(aaa, ' ', ' ', 5, atemps, nch)
    PRINT*, 'atemps est', atemps
    atemps = atemps(1:nch - 2)
    PRINT*, 'atemps', atemps
    read(atemps, *) imn
    dt = imn * 60
    PRINT*, 'le pas de temps dt', dt
    CALL getsch(aaa, ' ', ' ', 2, apasmax, nch)
    apasmax = apasmax(1:nch)
    read(apasmax, *) ipa
    pasmax = ipa
    PRINT*, 'pasmax est', pasmax
    CLOSE(97)
    !------------------------------------------------------------------
    ! *** on lit le pas de temps initial de la simulation ***
    !------------------------------------------------------------------
    in = itap
    !c                  pasprev=in
    !c                  time0=dt*(pasprev-1)
    pasprev = in - 1
    time0 = dt * pasprev

    close(98)

    WRITE(*, '(a)') 'OPEN ' // file_fordat
    open(99, FILE = file_fordat, FORM = 'UNFORMATTED', ACCESS = 'DIRECT', RECL = 8)
    icomp1 = nblvlm * 4
    IF (ts_fcg) icomp1 = icomp1 + 1
    IF (imp_fcg) icomp1 = icomp1 + nblvlm * 2
    IF (Turb_fcg) icomp1 = icomp1 + nblvlm * 2
    icompt = icomp1 * (in - 1)
    CALL rdgrads(99, icompt, nblvlm, z, ht_mes, hq_mes, hw_mes              &
            &, hu_mes, hv_mes, hthturb_mes, hqturb_mes              &
            &, ts_fcg, ts_subr, imp_fcg, Turb_fcg)
    PRINT *, 'get_uvd : rdgrads ->'
    PRINT *, tp_fcg

    ! following commented out because we have temperature already in ARM case
    !   (otherwise this is the potential temperature )
    !------------------------------------------------------------------------
    IF(Tp_fcg) THEN
      DO i = 1, nblvlm
        ht_mes(i) = ht_mes(i) * (hplaym(i) / 1000.)**rkappa
      enddo
    endif ! Tp_fcg
    PRINT*, 'ht_mes ', (ht_mes(i), i = 1, nblvlm)
    PRINT*, 'hq_mes ', (hq_mes(i), i = 1, nblvlm)
    PRINT*, 'hw_mes ', (hw_mes(i), i = 1, nblvlm)
    IF(imp_fcg) THEN
      PRINT*, 'hu_mes ', (hu_mes(i), i = 1, nblvlm)
      PRINT*, 'hv_mes ', (hv_mes(i), i = 1, nblvlm)
      PRINT*, 't', ts_subr
    endif ! imp_fcg
    IF(Turb_fcg) THEN
      PRINT*, 'hThTurb_mes ', (hThTurb_mes(i), i = 1, nblvlm)
      PRINT*, 'hqTurb ', (hqTurb_mes(i), i = 1, nblvlm)
    endif ! Turb_fcg
    !----------------------------------------------------------------------
    ! on a obtenu des champs initiaux sur les niveaux du meso_NH
    ! on interpole sur les niveaux du gcm(niveau pression bien sur!)
    !-----------------------------------------------------------------------
    DO k = 1, klev
      IF (JM(k) == 0) THEN
        !FKC bug? ne faut il pas convertir tsol en tendance ????
        !RT bug taken care of by removing the stuff
        htaft(k) = ht_mes(jm(k) + 1)
        hqaft(k) = hq_mes(jm(k) + 1)
        hwaft(k) = hw_mes(jm(k) + 1)
        IF(imp_fcg) THEN
          huaft(k) = hu_mes(jm(k) + 1)
          hvaft(k) = hv_mes(jm(k) + 1)
        endif ! imp_fcg
        IF(Turb_fcg) THEN
          hThTurbaft(k) = hThTurb_mes(jm(k) + 1)
          hqTurbaft(k) = hqTurb_mes(jm(k) + 1)
        endif ! Turb_fcg
      else ! JM(k) .EQ. 0
        htaft(k) = coef1(k) * ht_mes(jm(k)) + coef2(k) * ht_mes(jm(k) + 1)
        hqaft(k) = coef1(k) * hq_mes(jm(k)) + coef2(k) * hq_mes(jm(k) + 1)
        hwaft(k) = coef1(k) * hw_mes(jm(k)) + coef2(k) * hw_mes(jm(k) + 1)
        IF(imp_fcg) THEN
          huaft(k) = coef1(k) * hu_mes(jm(k)) + coef2(k) * hu_mes(jm(k) + 1)
          hvaft(k) = coef1(k) * hv_mes(jm(k)) + coef2(k) * hv_mes(jm(k) + 1)
        endif ! imp_fcg
        IF(Turb_fcg) THEN
          hThTurbaft(k) = coef1(k) * hThTurb_mes(jm(k))                        &
                  & + coef2(k) * hThTurb_mes(jm(k) + 1)
          hqTurbaft(k) = coef1(k) * hqTurb_mes(jm(k))                         &
                  & + coef2(k) * hqTurb_mes(jm(k) + 1)
        endif ! Turb_fcg
      endif ! JM(k) .EQ. 0
    enddo
    tsaft = ts_subr
    ! valeurs initiales des champs de convergence
    DO k = 1, klev
      ht(k) = htaft(k)
      hq(k) = hqaft(k)
      hw(k) = hwaft(k)
      IF(imp_fcg) THEN
        hu(k) = huaft(k)
        hv(k) = hvaft(k)
      endif ! imp_fcg
      IF(Turb_fcg) THEN
        hThTurb(k) = hThTurbaft(k)
        hqTurb(k) = hqTurbaft(k)
      endif ! Turb_fcg
    enddo
    ts_subr = tsaft
    close(99)
    close(98)

    !-------------------------------------------------------------------

    IF (Ts_fcg) Ts = Ts_subr
    return

    999     continue
    stop 'erreur lecture, file forcing.ctl'
  end

  SUBROUTINE advect_tvl(dtime, zt, zq, vu_f, vv_f, t_f, q_f                   &
          &, d_t_adv, d_q_adv)
    USE dimphy
USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
    IMPLICIT NONE


    !cccc      INCLUDE "dimphy.h"

    INTEGER k
    REAL dtime, fact, du, dv, cx, cy, alx, aly
    REAL zt(klev), zq(klev, 3)
    REAL vu_f(klev), vv_f(klev), t_f(klev), q_f(klev, 3)

    REAL d_t_adv(klev), d_q_adv(klev, 3)

    ! Velocity of moving cell
    data cx, cy /12., -2./

    ! Dimensions of moving cell
    data alx, aly /100000., 150000./

    DO k = 1, klev
      du = abs(vu_f(k) - cx) / alx
      dv = abs(vv_f(k) - cy) / aly
      fact = dtime * (du + dv - du * dv * dtime)
      d_t_adv(k) = fact * (t_f(k) - zt(k))
      d_q_adv(k, 1) = fact * (q_f(k, 1) - zq(k, 1))
      d_q_adv(k, 2) = fact * (q_f(k, 2) - zq(k, 2))
      d_q_adv(k, 3) = fact * (q_f(k, 3) - zq(k, 3))
    enddo

    RETURN
  end

  SUBROUTINE copie(klevgcm, playgcm, psolgcm, file_forctl)
    IMPLICIT NONE

    INTEGER k, klevgcm
    REAL playgcm(klevgcm) ! pression en milieu de couche du gcm
    REAL psolgcm
    CHARACTER*80 file_forctl

    klev = klevgcm

    !---------------------------------------------------------------------
    ! pression au milieu des couches du gcm dans la physiq
    ! (SB: remplace le CALL conv_lipress_gcm(playgcm) )
    !---------------------------------------------------------------------

    DO k = 1, klev
      play(k) = playgcm(k)
      PRINT*, 'la pression gcm est:', play(k)
    enddo

    !----------------------------------------------------------------------
    ! lecture du descripteur des donnees Meso-NH (forcing.ctl):
    !  -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
    !----------------------------------------------------------------------

    CALL mesolupbis(file_forctl)

    PRINT*, 'la valeur de nblvlm est:', nblvlm

    !----------------------------------------------------------------------
    ! etude de la correspondance entre les niveaux meso.NH et GCM;
    ! calcul des coefficients d interpolation coef1 et coef2
    !----------------------------------------------------------------------

    CALL corresbis(psolgcm)

    !---------------------------------------------------------
    ! TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
    !---------------------------------------------------------

    WRITE(*, *) ' '
    WRITE(*, *) 'TESTS com1_phys_gcss et com2_phys_gcss dans copie.F'
    WRITE(*, *) '--------------------------------------'
    WRITE(*, *) 'GCM: nb niveaux:', klev, ' et pression, coeffs:'
    DO k = 1, klev
      WRITE(*, *) play(k), coef1(k), coef2(k)
    enddo
    WRITE(*, *) 'MESO-NH: nb niveaux:', nblvlm, ' et pression:'
    DO k = 1, nblvlm
      WRITE(*, *) playm(k), hplaym(k)
    enddo
    WRITE(*, *) ' '

  END
  SUBROUTINE mesolupbis(file_forctl)
    !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

    ! Lecture descripteur des donnees MESO-NH (forcing.ctl):
    ! -------------------------------------------------------

    !     Cette subroutine lit dans le fichier de controle "essai.ctl"
    !     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
    !     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
    !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

    INTEGER i, lu, mlz, mlzh

    CHARACTER*80 file_forctl

    CHARACTER*4 a
    CHARACTER*80 aaa, anblvl
    INTEGER nch

    lu = 9
    open(lu, file = file_forctl, form = 'formatted')

    DO i = 1, 1000
      read(lu, 1000, end = 999) a
      IF (a == 'ZDEF') go to 100
    enddo

    100  backspace(lu)
    PRINT*, '  DESCRIPTION DES 2 MODELES : '
    PRINT*, ' '

    read(lu, 2000) aaa
    2000  format (a80)
    aaa = spaces(aaa, 1)
    CALL getsch(aaa, ' ', ' ', 2, anblvl, nch)
    read(anblvl, *) nblvlm

    PRINT*, 'nbre de niveaux de pression Meso-NH :', nblvlm
    PRINT*, ' '
    PRINT*, 'pression en Pa de chaque couche du meso-NH :'

    read(lu, *) (playm(mlz), mlz = 1, nblvlm)
    !      Si la pression est en HPa, la multiplier par 100
    IF (playm(1) < 10000.) THEN
      DO mlz = 1, nblvlm
        playm(mlz) = playm(mlz) * 100.
      enddo
    endif
    PRINT*, (playm(mlz), mlz = 1, nblvlm)

    1000 format (a4)

    PRINT*, ' '
    DO mlzh = 1, nblvlm
      hplaym(mlzh) = playm(mlzh) / 100.
    enddo

    PRINT*, 'pression en hPa de chaque couche du meso-NH: '
    PRINT*, (hplaym(mlzh), mlzh = 1, nblvlm)

    close (lu)
    return

    999  stop 'erreur lecture des niveaux pression des donnees'
  end

  SUBROUTINE rdgrads(itape, icount, nl, z, ht, hq, hw, hu, hv, hthtur, hqtur, &
          &  ts_fcg, ts, imp_fcg, Turb_fcg)
    IMPLICIT NONE
    INTEGER itape, icount, icomp, nl
    REAL z(nl), ht(nl), hq(nl), hw(nl), hu(nl), hv(nl)
    REAL hthtur(nl), hqtur(nl)
    REAL ts

    INTEGER k

    LOGICAL imp_fcg, ts_fcg, Turb_fcg

    icomp = icount

    DO k = 1, nl
      icomp = icomp + 1
      read(itape, rec = icomp)z(k)
      PRINT *, 'icomp,k,z(k) ', icomp, k, z(k)
    enddo
    DO k = 1, nl
      icomp = icomp + 1
      read(itape, rec = icomp)hT(k)
      PRINT*, hT(k), k
    enddo
    DO k = 1, nl
      icomp = icomp + 1
      read(itape, rec = icomp)hQ(k)
    enddo

    IF(turb_fcg) THEN
      DO k = 1, nl
        icomp = icomp + 1
        read(itape, rec = icomp)hThTur(k)
      enddo
      DO k = 1, nl
        icomp = icomp + 1
        read(itape, rec = icomp)hqTur(k)
      enddo
    endif
    PRINT *, ' apres lecture hthtur, hqtur'

    IF(imp_fcg) THEN
      DO k = 1, nl
        icomp = icomp + 1
        read(itape, rec = icomp)hu(k)
      enddo
      DO k = 1, nl
        icomp = icomp + 1
        read(itape, rec = icomp)hv(k)
      enddo

    endif

    DO k = 1, nl
      icomp = icomp + 1
      read(itape, rec = icomp)hw(k)
    enddo

    IF(ts_fcg) THEN
      icomp = icomp + 1
      read(itape, rec = icomp)ts
    endif

    PRINT *, ' rdgrads ->'

    RETURN
  END

  SUBROUTINE corresbis(psol)
    IMPLICIT NONE

    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    ! Cette subroutine calcule et affiche les valeurs des coefficients
    ! d interpolation qui serviront dans la formule d interpolation elle-
    ! meme.
    !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

    REAL psol
    REAL val
    INTEGER k, mlz

    DO k = 1, klev
      val = play(k)
      IF (val > playm(1)) THEN
        mlz = 0
        JM(1) = mlz
        coef1(1) = (playm(mlz + 1) - val) / (playm(mlz + 1) - psol)
        coef2(1) = (val - psol) / (playm(mlz + 1) - psol)
      ELSE IF (val > playm(nblvlm)) THEN
        DO mlz = 1, nblvlm
          IF (val <= playm(mlz).AND. val > playm(mlz + 1))THEN
            JM(k) = mlz
            coef1(k) = (playm(mlz + 1) - val) / (playm(mlz + 1) - playm(mlz))
            coef2(k) = (val - playm(mlz)) / (playm(mlz + 1) - playm(mlz))
          endif
        enddo
      else
        JM(k) = nblvlm - 1
        coef1(k) = 0.
        coef2(k) = 0.
      endif
    enddo

    !c      if (play(klev) .le. playm(nblvlm)) THEN
    !c         mlz=nblvlm-1
    !c         JM(klev)=mlz
    !c         coef1(klev)=(playm(mlz+1)-val)
    !c     *            /(playm(mlz+1)-playm(mlz))
    !c         coef2(klev)=(val-playm(mlz))
    !c     *            /(playm(mlz+1)-playm(mlz))
    !c      endif

    PRINT*, ' '
    PRINT*, '         INTERPOLATION  : '
    PRINT*, ' '
    PRINT*, 'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
    PRINT*, (JM(k), k = 1, klev)
    PRINT*, 'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
    PRINT*, (JM(k), k = 1, klev)
    PRINT*, ' '
    PRINT*, 'vals du premier coef d"interpolation pour les 9 niveaux: '
    PRINT*, (coef1(k), k = 1, klev)
    PRINT*, ' '
    PRINT*, 'valeurs du deuxieme coef d"interpolation pour les 9 niveaux:'
    PRINT*, (coef2(k), k = 1, klev)

    RETURN
  END
  SUBROUTINE GETSCH(STR, DEL, TRM, NTH, SST, NCH)
    !***************************************************************
    !*                                                             *
    !*                                                             *
    !* GETSCH                                                      *
    !*                                                             *
    !*                                                             *
    !* modified by :                                               *
    !***************************************************************
    !*   Return in SST the character string found between the NTH-1 and NTH
    !*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
    !*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
    !*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
    !*   NTH is greater than the number of delimiters in STR.
    IMPLICIT INTEGER (A-Z)
    CHARACTER STR*(*), DEL*1, TRM*1, SST*(*)
    NCH = -1
    SST = ' '
    IF(NTH>0) THEN
      IF(TRM==DEL) THEN
        LENGTH = LEN(STR)
      ELSE
        LENGTH = INDEX(STR, TRM) - 1
        IF(LENGTH<0) LENGTH = LEN(STR)
      ENDIF
      !*     Find beginning and end of the NTH DEL-limited substring in STR
      END = -1
      DO N = 1, NTH
        IF(END==LENGTH) RETURN
        BEG = END + 2
        END = BEG + INDEX(STR(BEG:LENGTH), DEL) - 2
        IF(END==BEG - 2) END = LENGTH
        !*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
      end DO
      NCH = END - BEG + 1
      IF(NCH>0) SST = STR(BEG:END)
    ENDIF
  END
  CHARACTER*(80) FUNCTION SPACES(STR, NSPACE)

    ! CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
    ! ORIG.  6/05/86 M.GOOSSENS/DD

    !-    The function value SPACES returns the character string STR with
    !-    leading blanks removed and each occurence of one or more blanks
    !-    replaced by NSPACE blanks inside the string STR

    CHARACTER*(80) STR
    INTEGER nspace, IBLANK, ISPACE, INONBL, LENSPA

    LENSPA = LEN(SPACES)
    SPACES = ' '
    IF (NSPACE<0) NSPACE = 0
    IBLANK = 1
    ISPACE = 1
    100 INONBL = INDEXC(STR(IBLANK:), ' ')
    IF (INONBL==0) THEN
      SPACES(ISPACE:) = STR(IBLANK:)
      RETURN
    ENDIF
    INONBL = INONBL + IBLANK - 1
    IBLANK = INDEX(STR(INONBL:), ' ')
    IF (IBLANK==0) THEN
      SPACES(ISPACE:) = STR(INONBL:)
      RETURN
    ENDIF
    IBLANK = IBLANK + INONBL - 1
    SPACES(ISPACE:) = STR(INONBL:IBLANK - 1)
    ISPACE = ISPACE + IBLANK - INONBL + NSPACE
    IF (ISPACE<=LENSPA)                         GO TO 100
  END
  INTEGER FUNCTION INDEXC(STR, SSTR)

    ! CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
    ! ORIG. 26/03/86 M.GOOSSENS/DD

    !-    Find the leftmost position where substring SSTR does not match
    !-    string STR scanning forward

    CHARACTER*(*) STR, SSTR
    INTEGER I, LENS, LENSS

    LENS = LEN(STR)
    LENSS = LEN(SSTR)

    DO I = 1, LENS - LENSS + 1
      IF (STR(I:I + LENSS - 1)/=SSTR) THEN
        INDEXC = I
        RETURN
      ENDIF
    END DO
    INDEXC = 0
  END
END MODULE lmdz_old_1dconv