Ignore:
Timestamp:
Apr 16, 2014, 7:16:58 AM (10 years ago)
Author:
fhourdin
Message:

Passage au format libre pour inclure les ficheirs du 1D dans lmdz1d.F90

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/1Dconv.h

    r2017 r2019  
    1         subroutine get_uvd(itap,dtime,file_forctl,file_fordat,
    2      :       ht,hq,hw,hu,hv,hthturb,hqturb,
    3      :       Ts,imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
    4 c
     1        subroutine get_uvd(itap,dtime,file_forctl,file_fordat,                  &
     2     &       ht,hq,hw,hu,hv,hthturb,hqturb,                                     &
     3     &       Ts,imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)                                 
     4!
    55        implicit none
    66 
    7 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    8 c cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
    9 c pouvoir calculer la convergence et le cisaillement dans la physiq
    10 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     7!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     8! cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
     9! pouvoir calculer la convergence et le cisaillement dans la physiq
     10!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    1111
    1212#include "YOMCST.h"
    … …  
    2929      COMMON/com2_phys_gcss/playm,hplaym,nblvlm
    3030
    31 c======================================================================
    32 c methode: on va chercher les donnees du mesoNH de meteo france, on y
    33 c          a acces a tout pas detemps grace a la routine rdgrads qui
    34 c          est une boucle lisant dans ces fichiers.
    35 c          Puis on interpole ces donnes sur les 11 niveaux du gcm et
    36 c          et sur les pas de temps de ce meme gcm
    37 c----------------------------------------------------------------------
    38 c input:
    39 c       pasmax     :nombre de pas de temps maximum du mesoNH
    40 c       dt         :pas de temps du meso_NH (en secondes)
    41 c----------------------------------------------------------------------
     31!======================================================================
     32! methode: on va chercher les donnees du mesoNH de meteo france, on y
     33!          a acces a tout pas detemps grace a la routine rdgrads qui
     34!          est une boucle lisant dans ces fichiers.
     35!          Puis on interpole ces donnes sur les 11 niveaux du gcm et
     36!          et sur les pas de temps de ce meme gcm
     37!----------------------------------------------------------------------
     38! input:
     39!       pasmax     :nombre de pas de temps maximum du mesoNH
     40!       dt         :pas de temps du meso_NH (en secondes)
     41!----------------------------------------------------------------------
    4242      integer pasmax,dt
    4343      save pasmax,dt
    44 c----------------------------------------------------------------------
    45 c arguments:
    46 c           itap   :compteur de la physique(le nombre de ces pas est
    47 c                   fixe dans la subroutine calcul_ini_gcm de interpo
    48 c                   -lation
    49 c           dtime  :pas detemps du gcm (en secondes)
    50 c           ht     :convergence horizontale de temperature(K/s)
    51 c           hq     :    "         "       d humidite (kg/kg/s)
    52 c           hw     :vitesse verticale moyenne (m/s**2)
    53 c           hu     :convergence horizontale d impulsion le long de x
    54 c                  (kg/(m^2 s^2)
    55 c           hv     : idem le long de y.
    56 c           Ts     : Temperature de surface (K)
    57 c           imp_fcg: var. logical .eq. T si forcage en impulsion
    58 c           ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier
    59 c           Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle
    60 c           Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier
    61 c----------------------------------------------------------------------
     44!----------------------------------------------------------------------
     45! arguments:
     46!           itap   :compteur de la physique(le nombre de ces pas est
     47!                   fixe dans la subroutine calcul_ini_gcm de interpo
     48!                   -lation
     49!           dtime  :pas detemps du gcm (en secondes)
     50!           ht     :convergence horizontale de temperature(K/s)
     51!           hq     :    "         "       d humidite (kg/kg/s)
     52!           hw     :vitesse verticale moyenne (m/s**2)
     53!           hu     :convergence horizontale d impulsion le long de x
     54!                  (kg/(m^2 s^2)
     55!           hv     : idem le long de y.
     56!           Ts     : Temperature de surface (K)
     57!           imp_fcg: var. logical .eq. T si forcage en impulsion
     58!           ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier
     59!           Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle
     60!           Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier
     61!----------------------------------------------------------------------
    6262        integer itap
    6363        real dtime
    … …  
    7474        logical Tp_fcg
    7575        logical Turb_fcg
    76 c----------------------------------------------------------------------
    77 c Variables internes de get_uvd (note : l interpolation temporelle
    78 c est faite entre les pas de temps before et after, sur les variables
    79 c definies sur la grille du SCM; on atteint exactement les valeurs Meso
    80 c aux milieux des pas de temps Meso)
    81 c     time0     :date initiale en secondes
    82 c     time      :temps associe a chaque pas du SCM
    83 c     pas       :numero du pas du meso_NH (on lit en pas : le premier pas
    84 c                 des donnees est duplique)
    85 c     pasprev   :numero du pas de lecture precedent
    86 c     htaft     :advection horizontale de temp. au pas de temps after
    87 c     hqaft     :    "         "      d humidite        "
    88 c     hwaft     :vitesse verticalle moyenne  au pas de temps after
    89 c     huaft,hvaft :advection horizontale d impulsion au pas de temps after
    90 c     tsaft     : surface temperature 'after time step'
    91 c     htbef     :idem htaft, mais pour le pas de temps before
    92 c     hqbef     :voir hqaft
    93 c     hwbef     :voir hwaft
    94 c     hubef,hvbef : idem huaft,hvaft, mais pour before
    95 c     tsbef     : surface temperature 'before time step'
    96 c----------------------------------------------------------------------
     76!----------------------------------------------------------------------
     77! Variables internes de get_uvd (note : l interpolation temporelle
     78! est faite entre les pas de temps before et after, sur les variables
     79! definies sur la grille du SCM; on atteint exactement les valeurs Meso
     80! aux milieux des pas de temps Meso)
     81!     time0     :date initiale en secondes
     82!     time      :temps associe a chaque pas du SCM
     83!     pas       :numero du pas du meso_NH (on lit en pas : le premier pas
     84!                 des donnees est duplique)
     85!     pasprev   :numero du pas de lecture precedent
     86!     htaft     :advection horizontale de temp. au pas de temps after
     87!     hqaft     :    "         "      d humidite        "
     88!     hwaft     :vitesse verticalle moyenne  au pas de temps after
     89!     huaft,hvaft :advection horizontale d impulsion au pas de temps after
     90!     tsaft     : surface temperature 'after time step'
     91!     htbef     :idem htaft, mais pour le pas de temps before
     92!     hqbef     :voir hqaft
     93!     hwbef     :voir hwaft
     94!     hubef,hvbef : idem huaft,hvaft, mais pour before
     95!     tsbef     : surface temperature 'before time step'
     96!----------------------------------------------------------------------
    9797        integer time0,pas,pasprev
    9898        save time0,pas,pasprev
    … …  
    106106        real Tsbef
    107107        save htbef,hqbef,hwbef,hubef,hvbef,hthturbbef,hqturbbef
    108 c
     108!
    109109        real timeaft,timebef
    110110        save timeaft,timebef
    111111        integer temps
    112112        character*4 string
    113 c----------------------------------------------------------------------
    114 c variables arguments de la subroutine rdgrads
    115 c---------------------------------------------------------------------
     113!----------------------------------------------------------------------
     114! variables arguments de la subroutine rdgrads
     115!---------------------------------------------------------------------
    116116        integer icompt,icomp1 !compteurs de rdgrads
    117117        real z(100)         ! altitude (grille Meso)
    … …  
    128128        real hqturb_mes(100) !tendance horizontale d humidite, due aux
    129129                              !flux turbulents
    130 c
    131 c---------------------------------------------------------------------
    132 c variable argument de la subroutine copie
    133 c---------------------------------------------------------------------
    134 c SB        real pplay(100)    !pression en milieu de couche du gcm
    135 c SB                            !argument de la physique
    136 c---------------------------------------------------------------------
    137 c variables destinees a la lecture du pas de temps du fichier de donnees
    138 c---------------------------------------------------------------------
     130!
     131!---------------------------------------------------------------------
     132! variable argument de la subroutine copie
     133!---------------------------------------------------------------------
     134! SB        real pplay(100)    !pression en milieu de couche du gcm
     135! SB                            !argument de la physique
     136!---------------------------------------------------------------------
     137! variables destinees a la lecture du pas de temps du fichier de donnees
     138!---------------------------------------------------------------------
    139139       character*80 aaa,atemps,spaces,apasmax
    140140       integer nch,imn,ipa
    141 c---------------------------------------------------------------------
    142 c  procedures appelees
     141!---------------------------------------------------------------------
     142!  procedures appelees
    143143        external rdgrads    !lire en iterant dans forcing.dat
    144 c---------------------------------------------------------------------
     144!---------------------------------------------------------------------
    145145               print*,'le pas itap est:',itap
    146 c*** on determine le pas du meso_NH correspondant au nouvel itap ***
    147 c*** pour aller chercher les champs dans rdgrads                 ***
    148 c
     146!*** on determine le pas du meso_NH correspondant au nouvel itap ***
     147!*** pour aller chercher les champs dans rdgrads                 ***
     148!
    149149        time=time0+itap*dtime
    150 cc        temps=int(time/dt+1)
    151 cc        pas=min(temps,pasmax)
     150!c        temps=int(time/dt+1)
     151!c        pas=min(temps,pasmax)
    152152        temps = 1 + int((dt + 2*time)/(2*dt))
    153153        pas=min(temps,pasmax-1)
    154154             print*,'le pas Meso est:',pas
    155 c
    156 c
    157 c===================================================================
    158 c
    159 c*** on remplit les champs before avec les champs after du pas   ***
    160 c*** precedent en format gcm                                     ***
     155!
     156!
     157!===================================================================
     158!
     159!*** on remplit les champs before avec les champs after du pas   ***
     160!*** precedent en format gcm                                     ***
    161161        if(pas.gt.pasprev)then
    162162          do i=1,klev
    … …  
    180180         print *, imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt
    181181                       print*,'le pas pas est:',pas
    182 c*** on va chercher les nouveaux champs after dans toga.dat     ***
    183 c*** champs en format meso_NH                                   ***
    184           open(99,FILE=file_fordat,FORM='UNFORMATTED',
    185      .             ACCESS='DIRECT',RECL=8)
    186           call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes
    187      .                  ,hu_mes,hv_mes,hthturb_mes,hqturb_mes
    188      .                  ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
    189 c
     182!*** on va chercher les nouveaux champs after dans toga.dat     ***
     183!*** champs en format meso_NH                                   ***
     184          open(99,FILE=file_fordat,FORM='UNFORMATTED',                        &
     185     &             ACCESS='DIRECT',RECL=8)
     186          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes                &
     187     &                  ,hu_mes,hv_mes,hthturb_mes,hqturb_mes                 &
     188     &                  ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
     189!
    190190
    191191               if(Tp_fcg) then
    192 c     (le forcage est donne en temperature potentielle)
     192!     (le forcage est donne en temperature potentielle)
    193193         do i = 1,nblvlm
    194194           ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
    … …  
    200200         enddo
    201201        endif  ! Turb_fcg
    202 c
     202!
    203203               print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
    204204               print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
    … …  
    213213                  endif
    214214          IF (ts_fcg) print*,'ts_subr', ts_subr
    215 c*** on interpole les champs meso_NH sur les niveaux de pression***
    216 c*** gcm . on obtient le nouveau champ after                    ***
     215!*** on interpole les champs meso_NH sur les niveaux de pression***
     216!*** gcm . on obtient le nouveau champ after                    ***
    217217            do k=1,klev
    218218             if (JM(k) .eq. 0) then
    … …  
    237237               endif ! imp_fcg
    238238               if(Turb_fcg) then
    239            hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))
    240      $               +coef2(k)*hThTurb_mes(jm(k)+1)
    241            hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))
    242      $               +coef2(k)*hqTurb_mes(jm(k)+1)
     239           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))                            &
     240     &               +coef2(k)*hThTurb_mes(jm(k)+1)
     241           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))                             &
     242     &               +coef2(k)*hqTurb_mes(jm(k)+1)
    243243               endif ! Turb_fcg
    244244             endif ! JM(k) .eq. 0
    … …  
    250250         endif ! pas.gt.pasprev    fin du bloc relatif au passage au pas
    251251                                  !de temps (meso) suivant
    252 c*** si on atteint le pas max des donnees experimentales ,on     ***
    253 c*** on conserve les derniers champs calcules                    ***
     252!*** si on atteint le pas max des donnees experimentales ,on     ***
     253!*** on conserve les derniers champs calcules                    ***
    254254      if(temps.ge.pasmax)then
    255255          do ll=1,klev
    … …  
    264264          ts_subr = tsaft
    265265      else ! temps.ge.pasmax
    266 c*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
    267 c** des pas de temps de 1h du meso_NH                            ***
     266!*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
     267!** des pas de temps de 1h du meso_NH                            ***
    268268         do j=1,klev
    269269         ht(j)=((timeaft-time)*htbef(j)+(time-timebef)*htaft(j))/dt
    … …  
    275275             endif ! imp_fcg
    276276             if(Turb_fcg) then
    277          hThTurb(j)=((timeaft-time)*hThTurbbef(j)
    278      $           +(time-timebef)*hThTurbaft(j))/dt
    279          hqTurb(j)= ((timeaft-time)*hqTurbbef(j)
    280      $           +(time-timebef)*hqTurbaft(j))/dt
     277         hThTurb(j)=((timeaft-time)*hThTurbbef(j)                           &
     278     &           +(time-timebef)*hThTurbaft(j))/dt
     279         hqTurb(j)= ((timeaft-time)*hqTurbbef(j)                            &
     280     &           +(time-timebef)*hqTurbaft(j))/dt
    281281             endif ! Turb_fcg
    282282         enddo
    283283         ts_subr = ((timeaft-time)*tsbef + (time-timebef)*tsaft)/dt
    284284       endif ! temps.ge.pasmax
    285 c
     285!
    286286        print *,' time,timebef,timeaft',time,timebef,timeaft
    287287        print *,' ht,htbef,htaft,hthturb,hthturbbef,hthturbaft'
    288288        do j= 1,klev
    289            print *, j,ht(j),htbef(j),htaft(j),
    290      $             hthturb(j),hthturbbef(j),hthturbaft(j)
     289           print *, j,ht(j),htbef(j),htaft(j),                              &
     290     &             hthturb(j),hthturbbef(j),hthturbaft(j)
    291291        enddo
    292292        print *,' hq,hqbef,hqaft,hqturb,hqturbbef,hqturbaft'
    293293        do j= 1,klev
    294            print *, j,hq(j),hqbef(j),hqaft(j),
    295      $             hqturb(j),hqturbbef(j),hqturbaft(j)
     294           print *, j,hq(j),hqbef(j),hqaft(j),                              &
     295     &             hqturb(j),hqturbbef(j),hqturbaft(j)
    296296        enddo
    297 c
    298 c-------------------------------------------------------------------
    299 c
     297!
     298!-------------------------------------------------------------------
     299!
    300300         IF (Ts_fcg) Ts = Ts_subr
    301301         return
    302 c
    303 c-----------------------------------------------------------------------
    304 c on sort les champs de "convergence" pour l instant initial 'in'
    305 c ceci se passe au pas temps itap=0 de la physique
    306 c-----------------------------------------------------------------------
    307         entry get_uvd2(itap,dtime,file_forctl,file_fordat,
    308      .           ht,hq,hw,hu,hv,hThTurb,hqTurb,ts,
    309      .           imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
     302!
     303!-----------------------------------------------------------------------
     304! on sort les champs de "convergence" pour l instant initial 'in'
     305! ceci se passe au pas temps itap=0 de la physique
     306!-----------------------------------------------------------------------
     307        entry get_uvd2(itap,dtime,file_forctl,file_fordat,                  &
     308     &           ht,hq,hw,hu,hv,hThTurb,hqTurb,ts,                          &
     309     &           imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
    310310             print*,'le pas itap est:',itap
    311 c
    312 c===================================================================
    313 c
     311!
     312!===================================================================
     313!
    314314       write(*,'(a)') 'OPEN '//file_forctl
    315315       open(97,FILE=file_forctl,FORM='FORMATTED')
    316 c
    317 c------------------
     316!
     317!------------------
    318318      do i=1,1000
    319319      read(97,1000,end=999) string
    … …  
    322322      enddo
    323323 50   backspace(97)
    324 c-------------------------------------------------------------------
    325 c   *** on lit le pas de temps dans le fichier de donnees ***
    326 c   *** "forcing.ctl" et pasmax                           ***
    327 c-------------------------------------------------------------------
     324!-------------------------------------------------------------------
     325!   *** on lit le pas de temps dans le fichier de donnees ***
     326!   *** "forcing.ctl" et pasmax                           ***
     327!-------------------------------------------------------------------
    328328      read(97,2000) aaa
    329329 2000  format (a80)
    … …  
    344344         print*,'pasmax est',pasmax
    345345      CLOSE(97)
    346 c------------------------------------------------------------------
    347 c *** on lit le pas de temps initial de la simulation ***
    348 c------------------------------------------------------------------
     346!------------------------------------------------------------------
     347! *** on lit le pas de temps initial de la simulation ***
     348!------------------------------------------------------------------
    349349                  in=itap
    350 cc                  pasprev=in
    351 cc                  time0=dt*(pasprev-1)
     350!c                  pasprev=in
     351!c                  time0=dt*(pasprev-1)
    352352                  pasprev=in-1
    353353                  time0=dt*pasprev
    354 C
     354!
    355355          close(98)
    356 c
     356!
    357357      write(*,'(a)') 'OPEN '//file_fordat
    358       open(99,FILE=file_fordat,FORM='UNFORMATTED',
    359      .          ACCESS='DIRECT',RECL=8)
     358      open(99,FILE=file_fordat,FORM='UNFORMATTED',                          &
     359     &          ACCESS='DIRECT',RECL=8)
    360360          icomp1 = nblvlm*4
    361361          IF (ts_fcg) icomp1 = icomp1 + 1
    … …  
    363363          IF (Turb_fcg) icomp1 = icomp1 + nblvlm*2
    364364          icompt = icomp1*(in-1)
    365           call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes
    366      .                   ,hu_mes,hv_mes,hthturb_mes,hqturb_mes
    367      .                   ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
     365          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes              &
     366     &                   ,hu_mes,hv_mes,hthturb_mes,hqturb_mes              &
     367     &                   ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
    368368          print *, 'get_uvd : rdgrads ->'
    369369          print *, tp_fcg
    370 c
    371 c following commented out because we have temperature already in ARM case
    372 c   (otherwise this is the potential temperature )
    373 c------------------------------------------------------------------------
     370!
     371! following commented out because we have temperature already in ARM case
     372!   (otherwise this is the potential temperature )
     373!------------------------------------------------------------------------
    374374               if(Tp_fcg) then
    375375          do i = 1,nblvlm
    … …  
    389389                 print*,'hqTurb ',     (hqTurb_mes(i),i=1,nblvlm)
    390390              endif ! Turb_fcg
    391 c----------------------------------------------------------------------
    392 c on a obtenu des champs initiaux sur les niveaux du meso_NH
    393 c on interpole sur les niveaux du gcm(niveau pression bien sur!)
    394 c-----------------------------------------------------------------------
     391!----------------------------------------------------------------------
     392! on a obtenu des champs initiaux sur les niveaux du meso_NH
     393! on interpole sur les niveaux du gcm(niveau pression bien sur!)
     394!-----------------------------------------------------------------------
    395395            do k=1,klev
    396396             if (JM(k) .eq. 0) then
    397 cFKC bug? ne faut il pas convertir tsol en tendance ????
    398 cRT bug taken care of by removing the stuff
     397!FKC bug? ne faut il pas convertir tsol en tendance ????
     398!RT bug taken care of by removing the stuff
    399399           htaft(k)=              ht_mes(jm(k)+1)
    400400           hqaft(k)=              hq_mes(jm(k)+1)
    … …  
    417417               endif ! imp_fcg
    418418               if(Turb_fcg) then
    419            hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))
    420      $                  +coef2(k)*hThTurb_mes(jm(k)+1)
    421            hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))
    422      $                  +coef2(k)*hqTurb_mes(jm(k)+1)
     419           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))                        &
     420     &                  +coef2(k)*hThTurb_mes(jm(k)+1)
     421           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))                         &
     422     &                  +coef2(k)*hqTurb_mes(jm(k)+1)
    423423               endif ! Turb_fcg
    424424             endif ! JM(k) .eq. 0
    425425            enddo
    426426            tsaft = ts_subr
    427 c valeurs initiales des champs de convergence
     427! valeurs initiales des champs de convergence
    428428          do k=1,klev
    429429             ht(k)=htaft(k)
    … …  
    442442          close(99)
    443443          close(98)
    444 c
    445 c-------------------------------------------------------------------
    446 c
    447 c
     444!
     445!-------------------------------------------------------------------
     446!
     447!
    448448 100      IF (Ts_fcg) Ts = Ts_subr
    449449        return
    450 c
     450!
    451451999     continue
    452452        stop 'erreur lecture, file forcing.ctl'
    453453        end
    454454
    455       SUBROUTINE advect_tvl(dtime,zt,zq,vu_f,vv_f,t_f,q_f
    456      :                     ,d_t_adv,d_q_adv)
     455      SUBROUTINE advect_tvl(dtime,zt,zq,vu_f,vv_f,t_f,q_f                   &
     456     &                     ,d_t_adv,d_q_adv)
    457457      use dimphy
    458458      implicit none
    459459
    460460#include "dimensions.h"
    461 ccccc#include "dimphy.h"
     461!cccc#include "dimphy.h"
    462462
    463463      integer k
    464464      real dtime, fact, du, dv, cx, cy, alx, aly
    465465      real zt(klev), zq(klev,3)
    466      :   , vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
     466      real vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
    467467
    468468      real d_t_adv(klev), d_q_adv(klev,3)
    469469
    470 c Velocity of moving cell
     470! Velocity of moving cell
    471471      data cx,cy /12., -2./
    472472
    473 c Dimensions of moving cell
    474       data alx,aly /100 000.,150 000./
     473! Dimensions of moving cell
     474      data alx,aly /100000.,150000./
    475475
    476476      do k = 1, klev
    … …  
    490490      implicit none
    491491
    492 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    493 c cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h
    494 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     492!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     493! cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h
     494!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    495495
    496496      INTEGER klev !nombre de niveau de pression du GCM
    … …  
    514514      klev = klevgcm
    515515
    516 c---------------------------------------------------------------------
    517 c pression au milieu des couches du gcm dans la physiq
    518 c (SB: remplace le call conv_lipress_gcm(playgcm) )
    519 c---------------------------------------------------------------------
     516!---------------------------------------------------------------------
     517! pression au milieu des couches du gcm dans la physiq
     518! (SB: remplace le call conv_lipress_gcm(playgcm) )
     519!---------------------------------------------------------------------
    520520
    521521       do k = 1, klev
    … …  
    524524       enddo
    525525
    526 c----------------------------------------------------------------------
    527 c lecture du descripteur des donnees Meso-NH (forcing.ctl):
    528 c  -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
    529 c (on remplit le COMMON com2_phys_gcss)
    530 c----------------------------------------------------------------------
     526!----------------------------------------------------------------------
     527! lecture du descripteur des donnees Meso-NH (forcing.ctl):
     528!  -> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
     529! (on remplit le COMMON com2_phys_gcss)
     530!----------------------------------------------------------------------
    531531
    532532      call mesolupbis(file_forctl)
    … …  
    534534      print*,'la valeur de nblvlm est:',nblvlm
    535535
    536 c----------------------------------------------------------------------
    537 c etude de la correspondance entre les niveaux meso.NH et GCM;
    538 c calcul des coefficients d interpolation coef1 et coef2
    539 c (on remplit le COMMON com1_phys_gcss)
    540 c----------------------------------------------------------------------
     536!----------------------------------------------------------------------
     537! etude de la correspondance entre les niveaux meso.NH et GCM;
     538! calcul des coefficients d interpolation coef1 et coef2
     539! (on remplit le COMMON com1_phys_gcss)
     540!----------------------------------------------------------------------
    541541
    542542      call corresbis(psolgcm)
    543543
    544 c---------------------------------------------------------
    545 c TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
    546 c---------------------------------------------------------
     544!---------------------------------------------------------
     545! TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
     546!---------------------------------------------------------
    547547 
    548548      write(*,*) ' '
    … …  
    562562      SUBROUTINE mesolupbis(file_forctl)
    563563      implicit none
    564 c
    565 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    566 c
    567 c Lecture descripteur des donnees MESO-NH (forcing.ctl):
    568 c -------------------------------------------------------
    569 c
    570 c     Cette subroutine lit dans le fichier de controle "essai.ctl"
    571 c     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
    572 c     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
    573 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    574 c
     564!
     565!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     566!
     567! Lecture descripteur des donnees MESO-NH (forcing.ctl):
     568! -------------------------------------------------------
     569!
     570!     Cette subroutine lit dans le fichier de controle "essai.ctl"
     571!     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
     572!     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
     573!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     574!
    575575      INTEGER nblvlm !nombre de niveau de pression du mesoNH
    576576      REAL playm(100)  !pression en Pa milieu de chaque couche Meso-NH
    … …  
    588588      lu=9
    589589      open(lu,file=file_forctl,form='formatted')
    590 c
     590!
    591591      do i=1,1000
    592592      read(lu,1000,end=999) a
    593593      if (a .eq. 'ZDEF') go to 100
    594594      enddo
    595 c
     595!
    596596 100  backspace(lu)
    597597      print*,'  DESCRIPTION DES 2 MODELES : '
    598598      print*,' '
    599 c
     599!
    600600      read(lu,2000) aaa
    601601 2000  format (a80)
    … …  
    604604         read(anblvl,*) nblvlm
    605605
    606 c
     606!
    607607      print*,'nbre de niveaux de pression Meso-NH :',nblvlm
    608608      print*,' '
    609609      print*,'pression en Pa de chaque couche du meso-NH :'
    610 c
     610!
    611611      read(lu,*) (playm(mlz),mlz=1,nblvlm)
    612 c      Si la pression est en HPa, la multiplier par 100
     612!      Si la pression est en HPa, la multiplier par 100
    613613      if (playm(1) .lt. 10000.) then
    614614        do mlz = 1,nblvlm
    … …  
    617617      endif
    618618      print*,(playm(mlz),mlz=1,nblvlm)
    619 c
     619!
    620620 1000 format (a4)
    621621 1001 format(5x,i2)
    622 c
     622!
    623623      print*,' '
    624624      do mlzh=1,nblvlm
    625625      hplaym(mlzh)=playm(mlzh)/100.
    626626      enddo
    627 c
     627!
    628628      print*,'pression en hPa de chaque couche du meso-NH: '
    629629      print*,(hplaym(mlzh),mlzh=1,nblvlm)
    630 c
     630!
    631631      close (lu)
    632632      return
    633 c
     633!
    634634 999  stop 'erreur lecture des niveaux pression des donnees'
    635635      end
    636636
    637       SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw,hu,hv,hthtur,hqtur,
    638      :  ts_fcg,ts,imp_fcg,Turb_fcg)
     637      SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw,hu,hv,hthtur,hqtur,     &
     638     &  ts_fcg,ts,imp_fcg,Turb_fcg)
    639639      IMPLICIT none
    640640      INTEGER itape,icount,icomp, nl
    … …  
    642642      real hthtur(nl),hqtur(nl)
    643643      real ts
    644 c
     644!
    645645      INTEGER k
    646 c
     646!
    647647      LOGICAL imp_fcg,ts_fcg,Turb_fcg
    648 c
     648!
    649649      icomp = icount
    650 c
    651 c
     650!
     651!
    652652         do k=1,nl
    653653            icomp=icomp+1
    … …  
    664664            read(itape,rec=icomp)hQ(k)
    665665         enddo
    666 c
     666!
    667667             if(turb_fcg) then
    668668         do k=1,nl
    … …  
    676676             endif
    677677         print *,' apres lecture hthtur, hqtur'
    678 c
     678!
    679679          if(imp_fcg) then
    680680
    … …  
    689689
    690690          endif
    691 c
     691!
    692692         do k=1,nl
    693693            icomp=icomp+1
    694694            read(itape,rec=icomp)hw(k)
    695695         enddo
    696 c
     696!
    697697              if(ts_fcg) then
    698698         icomp=icomp+1
    699699         read(itape,rec=icomp)ts
    700700              endif
    701 c
     701!
    702702      print *,' rdgrads ->'
    703703
    … …  
    708708      implicit none
    709709
    710 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    711 c Cette subroutine calcule et affiche les valeurs des coefficients
    712 c d interpolation qui serviront dans la formule d interpolation elle-
    713 c meme.
    714 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     710!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     711! Cette subroutine calcule et affiche les valeurs des coefficients
     712! d interpolation qui serviront dans la formule d interpolation elle-
     713! meme.
     714!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    715715
    716716      INTEGER klev    !nombre de niveau de pression du GCM
    … …  
    737737          mlz = 0
    738738          JM(1) = mlz
    739           coef1(1)=(playm(mlz+1)-val)
    740      *             /(playm(mlz+1)-psol)
    741           coef2(1)=(val-psol)
    742      *             /(playm(mlz+1)-psol)
     739          coef1(1)=(playm(mlz+1)-val)/(playm(mlz+1)-psol)
     740          coef2(1)=(val-psol)/(playm(mlz+1)-psol)
    743741       else if (val .gt. playm(nblvlm)) then
    744742         do mlz=1,nblvlm
    745           if (     val .le. playm(mlz)
    746      *       .and. val .gt. playm(mlz+1))then
     743          if (     val .le. playm(mlz).and. val .gt. playm(mlz+1))then
    747744           JM(k)=mlz
    748            coef1(k)=(playm(mlz+1)-val)
    749      *              /(playm(mlz+1)-playm(mlz))
    750            coef2(k)=(val-playm(mlz))
    751      *              /(playm(mlz+1)-playm(mlz))
     745           coef1(k)=(playm(mlz+1)-val)/(playm(mlz+1)-playm(mlz))
     746           coef2(k)=(val-playm(mlz))/(playm(mlz+1)-playm(mlz))
    752747          endif
    753748         enddo
    … …  
    758753       endif
    759754      enddo
    760 c
    761 cc      if (play(klev) .le. playm(nblvlm)) then
    762 cc         mlz=nblvlm-1
    763 cc         JM(klev)=mlz
    764 cc         coef1(klev)=(playm(mlz+1)-val)
    765 cc     *            /(playm(mlz+1)-playm(mlz))
    766 cc         coef2(klev)=(val-playm(mlz))
    767 cc     *            /(playm(mlz+1)-playm(mlz))
    768 cc      endif
    769 c
     755!
     756!c      if (play(klev) .le. playm(nblvlm)) then
     757!c         mlz=nblvlm-1
     758!c         JM(klev)=mlz
     759!c         coef1(klev)=(playm(mlz+1)-val)
     760!c     *            /(playm(mlz+1)-playm(mlz))
     761!c         coef2(klev)=(val-playm(mlz))
     762!c     *            /(playm(mlz+1)-playm(mlz))
     763!c      endif
     764!
    770765      print*,' '
    771766      print*,'         INTERPOLATION  : '
    … …  
    776771      print*,(JM(k),k=1,klev)
    777772      print*,' '
    778       print*,'valeurs du premier coef d"interpolation pour les 9 niveaux
    779      *: '
     773      print*,'vals du premier coef d"interpolation pour les 9 niveaux: '
    780774      print*,(coef1(k),k=1,klev)
    781775      print*,' '
    782       print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveau
    783      *x: '
     776      print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveaux:'
    784777      print*,(coef2(k),k=1,klev)
    785 c
     778!
    786779      return
    787780      end
    788781      SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH)
    789 C***************************************************************
    790 C*                                                             *
    791 C*                                                             *
    792 C* GETSCH                                                      *
    793 C*                                                             *
    794 C*                                                             *
    795 C* modified by :                                               *
    796 C***************************************************************
    797 C*   Return in SST the character string found between the NTH-1 and NTH
    798 C*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
    799 C*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
    800 C*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
    801 C*   NTH is greater than the number of delimiters in STR.
     782!***************************************************************
     783!*                                                             *
     784!*                                                             *
     785!* GETSCH                                                      *
     786!*                                                             *
     787!*                                                             *
     788!* modified by :                                               *
     789!***************************************************************
     790!*   Return in SST the character string found between the NTH-1 and NTH
     791!*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
     792!*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
     793!*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
     794!*   NTH is greater than the number of delimiters in STR.
    802795      IMPLICIT INTEGER (A-Z)
    803796      CHARACTER STR*(*),DEL*1,TRM*1,SST*(*)
    … …  
    811804          IF(LENGTH.LT.0) LENGTH=LEN(STR)
    812805        ENDIF
    813 C*     Find beginning and end of the NTH DEL-limited substring in STR
     806!*     Find beginning and end of the NTH DEL-limited substring in STR
    814807        END=-1
    815808        DO 1,N=1,NTH
    … …  
    818811        END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2
    819812        IF(END.EQ.BEG-2) END=LENGTH
    820 C*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
     813!*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
    821814    1   CONTINUE
    822815        NCH=END-BEG+1
    … …  
    825818      END
    826819      CHARACTER*(*) FUNCTION SPACES(STR,NSPACE)
    827 C
    828 C CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
    829 C ORIG.  6/05/86 M.GOOSSENS/DD
    830 C
    831 C-    The function value SPACES returns the character string STR with
    832 C-    leading blanks removed and each occurence of one or more blanks
    833 C-    replaced by NSPACE blanks inside the string STR
    834 C
     820!
     821! CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
     822! ORIG.  6/05/86 M.GOOSSENS/DD
     823!
     824!-    The function value SPACES returns the character string STR with
     825!-    leading blanks removed and each occurence of one or more blanks
     826!-    replaced by NSPACE blanks inside the string STR
     827!
    835828      CHARACTER*(*) STR
    836 C
     829!
    837830      LENSPA = LEN(SPACES)
    838831      SPACES = ' '
    … …  
    857850  999 END
    858851      FUNCTION INDEXC(STR,SSTR)
    859 C
    860 C CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
    861 C ORIG. 26/03/86 M.GOOSSENS/DD
    862 C
    863 C-    Find the leftmost position where substring SSTR does not match
    864 C-    string STR scanning forward
    865 C
     852!
     853! CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
     854! ORIG. 26/03/86 M.GOOSSENS/DD
     855!
     856!-    Find the leftmost position where substring SSTR does not match
     857!-    string STR scanning forward
     858!
    866859      CHARACTER*(*) STR,SSTR
    867 C
     860!
    868861      LENS   = LEN(STR)
    869862      LENSS  = LEN(SSTR)
    870 C
     863!
    871864      DO 10 I=1,LENS-LENSS+1
    872865          IF (STR(I:I+LENSS-1).NE.SSTR) THEN
    … …  
    876869   10 CONTINUE
    877870      INDEXC = 0
    878 C
     871!
    879872  999 END
Note: See TracChangeset for help on using the changeset viewer.