source: LMDZ5/branches/LF-private/libf/phy1d/1Dconv.h @ 5405

Last change on this file since 5405 was 1607, checked in by lguez, 13 years ago

Import 1D files. Files added to directory "phy1d" were in directories:

lmdz1d_source_20111207/phy1d_source
lmdz1d_source_20111207/phy1d_source_upd

extracted from:

http://www.lmd.jussieu.fr/~jyg/lmdz1d_source_20111207.tar.gz

  • Property svn:executable set to *
File size: 31.4 KB
Line 
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)
4c
5        implicit none
6 
7cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
8c cette routine permet d obtenir u_convg,v_convg,ht,hq et ainsi de
9c pouvoir calculer la convergence et le cisaillement dans la physiq
10ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
11
12#include "YOMCST.h"
13
14      INTEGER klev
15      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
16      INTEGER JM(100) !pression en Pa au milieu de chaque couche GCM
17      REAL coef1(100) !coefficient d interpolation
18      REAL coef2(100) !coefficient d interpolation
19
20      INTEGER nblvlm !nombre de niveau de pression du mesoNH
21      REAL playm(100)  !pression en Pa au milieu de chaque couche Meso-NH
22      REAL hplaym(100) !pression en hPa milieux des couches Meso-NH
23
24      integer i,j,k,ii,ll,in
25      REAL tsol,qsol
26
27      CHARACTER*80 file_forctl,file_fordat
28
29      COMMON/com1_phys_gcss/klev,play,JM,coef1,coef2
30      COMMON/com2_phys_gcss/nblvlm,playm,hplaym
31
32c======================================================================
33c methode: on va chercher les donnees du mesoNH de meteo france, on y
34c          a acces a tout pas detemps grace a la routine rdgrads qui
35c          est une boucle lisant dans ces fichiers.
36c          Puis on interpole ces donnes sur les 11 niveaux du gcm et
37c          et sur les pas de temps de ce meme gcm
38c----------------------------------------------------------------------
39c input:
40c       pasmax     :nombre de pas de temps maximum du mesoNH
41c       dt         :pas de temps du meso_NH (en secondes)
42c----------------------------------------------------------------------
43      integer pasmax,dt
44      save pasmax,dt
45c----------------------------------------------------------------------
46c arguments:
47c           itap   :compteur de la physique(le nombre de ces pas est
48c                   fixe dans la subroutine calcul_ini_gcm de interpo
49c                   -lation
50c           dtime  :pas detemps du gcm (en secondes)
51c           ht     :convergence horizontale de temperature(K/s)
52c           hq     :    "         "       d humidite (kg/kg/s)
53c           hw     :vitesse verticale moyenne (m/s**2)
54c           hu     :convergence horizontale d impulsion le long de x
55c                  (kg/(m^2 s^2)
56c           hv     : idem le long de y.
57c           Ts     : Temperature de surface (K)
58c           imp_fcg: var. logical .eq. T si forcage en impulsion
59c           ts_fcg: var. logical .eq. T si forcage en Ts present dans fichier
60c           Tp_fcg: var. logical .eq. T si forcage donne en Temp potentielle
61c           Turb_fcg: var. logical .eq. T si forcage turbulent present dans fichier
62c----------------------------------------------------------------------
63        integer itap
64        real dtime
65        real ht(100)
66        real hq(100)
67        real hu(100)
68        real hv(100)
69        real hw(100)
70        real hthturb(100)
71        real hqturb(100)
72        real Ts, Ts_subr
73        logical imp_fcg
74        logical ts_fcg
75        logical Tp_fcg
76        logical Turb_fcg
77c----------------------------------------------------------------------
78c Variables internes de get_uvd (note : l interpolation temporelle
79c est faite entre les pas de temps before et after, sur les variables
80c definies sur la grille du SCM; on atteint exactement les valeurs Meso
81c aux milieux des pas de temps Meso)
82c     time0     :date initiale en secondes
83c     time      :temps associe a chaque pas du SCM
84c     pas       :numero du pas du meso_NH (on lit en pas : le premier pas
85c                 des donnees est duplique)
86c     pasprev   :numero du pas de lecture precedent
87c     htaft     :advection horizontale de temp. au pas de temps after
88c     hqaft     :    "         "      d humidite        "
89c     hwaft     :vitesse verticalle moyenne  au pas de temps after
90c     huaft,hvaft :advection horizontale d impulsion au pas de temps after
91c     tsaft     : surface temperature 'after time step'
92c     htbef     :idem htaft, mais pour le pas de temps before
93c     hqbef     :voir hqaft
94c     hwbef     :voir hwaft
95c     hubef,hvbef : idem huaft,hvaft, mais pour before
96c     tsbef     : surface temperature 'before time step'
97c----------------------------------------------------------------------
98        integer time0,pas,pasprev
99        save time0,pas,pasprev
100        real time
101        real htaft(100),hqaft(100),hwaft(100),huaft(100),hvaft(100)
102        real hthturbaft(100),hqturbaft(100)
103        real Tsaft
104        save htaft,hqaft,hwaft,huaft,hvaft,hthturbaft,hqturbaft
105        real htbef(100),hqbef(100),hwbef(100),hubef(100),hvbef(100)
106        real hthturbbef(100),hqturbbef(100)
107        real Tsbef
108        save htbef,hqbef,hwbef,hubef,hvbef,hthturbbef,hqturbbef
109c
110        real timeaft,timebef
111        save timeaft,timebef
112        integer temps
113        character*4 string
114c----------------------------------------------------------------------
115c variables arguments de la subroutine rdgrads
116c---------------------------------------------------------------------
117        integer icompt,icomp1 !compteurs de rdgrads
118        real z(100)         ! altitude (grille Meso)
119        real ht_mes(100)    !convergence horizontale de temperature
120                            !-(grille Meso)
121        real hq_mes(100)    !convergence horizontale d humidite
122                            !(grille Meso)
123        real hw_mes(100)    !vitesse verticale moyenne
124                            !(grille Meso)
125        real hu_mes(100),hv_mes(100)    !convergence horizontale d impulsion
126                                        !(grille Meso)
127        real hthturb_mes(100) !tendance horizontale de T_pot, due aux
128                              !flux turbulents
129        real hqturb_mes(100) !tendance horizontale d humidite, due aux
130                              !flux turbulents
131c
132c---------------------------------------------------------------------
133c variable argument de la subroutine copie
134c---------------------------------------------------------------------
135c SB        real pplay(100)    !pression en milieu de couche du gcm
136c SB                            !argument de la physique
137c---------------------------------------------------------------------
138c variables destinees a la lecture du pas de temps du fichier de donnees
139c---------------------------------------------------------------------
140       character*80 aaa,atemps,spaces,apasmax
141       integer nch,imn,ipa
142c---------------------------------------------------------------------
143c  procedures appelees
144        external rdgrads    !lire en iterant dans forcing.dat
145c---------------------------------------------------------------------
146               print*,'le pas itap est:',itap
147c*** on determine le pas du meso_NH correspondant au nouvel itap ***
148c*** pour aller chercher les champs dans rdgrads                 ***
149c
150        time=time0+itap*dtime
151cc        temps=int(time/dt+1)
152cc        pas=min(temps,pasmax)
153        temps = 1 + int((dt + 2*time)/(2*dt))
154        pas=min(temps,pasmax-1)
155             print*,'le pas Meso est:',pas
156c
157c
158c===================================================================
159c
160c*** on remplit les champs before avec les champs after du pas   ***
161c*** precedent en format gcm                                     ***
162        if(pas.gt.pasprev)then
163          do i=1,klev
164             htbef(i)=htaft(i)
165             hqbef(i)=hqaft(i)
166             hwbef(i)=hwaft(i)
167             hubef(i)=huaft(i)
168             hvbef(i)=hvaft(i)
169             hThTurbbef(i)=hThTurbaft(i)
170             hqTurbbef(i)=hqTurbaft(i)
171          enddo
172               tsbef = tsaft
173          timebef=pasprev*dt
174          timeaft=timebef+dt
175          icomp1 = nblvlm*4
176          IF (ts_fcg) icomp1 = icomp1 + 1
177          IF (imp_fcg) icomp1 = icomp1 + nblvlm*2
178          IF (Turb_fcg) icomp1 = icomp1 + nblvlm*2
179          icompt = icomp1*pas
180         print *, 'imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt'
181         print *, imp_fcg,ts_fcg,Turb_fcg,pas,nblvlm,icompt
182                       print*,'le pas pas est:',pas
183c*** on va chercher les nouveaux champs after dans toga.dat     ***
184c*** champs en format meso_NH                                   ***
185          open(99,FILE=file_fordat,FORM='UNFORMATTED',
186     .             ACCESS='DIRECT',RECL=8)
187          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes
188     .                  ,hu_mes,hv_mes,hthturb_mes,hqturb_mes
189     .                  ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
190c
191
192               if(Tp_fcg) then
193c     (le forcage est donne en temperature potentielle)
194         do i = 1,nblvlm
195           ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
196         enddo
197               endif ! Tp_fcg
198        if(Turb_fcg) then
199         do i = 1,nblvlm
200           hThTurb_mes(i) = hThTurb_mes(i)*(hplaym(i)/1000.)**rkappa
201         enddo
202        endif  ! Turb_fcg
203c
204               print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
205               print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
206               print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
207                  if(imp_fcg) then
208               print*,'hu_mes ',(hu_mes(i),i=1,nblvlm)
209               print*,'hv_mes ',(hv_mes(i),i=1,nblvlm)
210                  endif
211                  if(Turb_fcg) then
212               print*,'hThTurb_mes ',(hThTurb_mes(i),i=1,nblvlm)
213               print*,'hqTurb_mes ',(hqTurb_mes(i),i=1,nblvlm)
214                  endif
215               IF (ts_fcg) print*,'ts_subr', ts_subr
216c*** on interpole les champs meso_NH sur les niveaux de pression***
217c*** gcm . on obtient le nouveau champ after                    ***
218            do k=1,klev
219             if (JM(k) .eq. 0) then
220         htaft(k)=              ht_mes(jm(k)+1)
221         hqaft(k)=              hq_mes(jm(k)+1)
222         hwaft(k)=              hw_mes(jm(k)+1)
223               if(imp_fcg) then
224           huaft(k)=              hu_mes(jm(k)+1)
225           hvaft(k)=              hv_mes(jm(k)+1)
226               endif ! imp_fcg
227               if(Turb_fcg) then
228           hThTurbaft(k)=         hThTurb_mes(jm(k)+1)
229           hqTurbaft(k)=          hqTurb_mes(jm(k)+1)
230               endif ! Turb_fcg
231             else ! JM(k) .eq. 0
232           htaft(k)=coef1(k)*ht_mes(jm(k))+coef2(k)*ht_mes(jm(k)+1)
233           hqaft(k)=coef1(k)*hq_mes(jm(k))+coef2(k)*hq_mes(jm(k)+1)
234           hwaft(k)=coef1(k)*hw_mes(jm(k))+coef2(k)*hw_mes(jm(k)+1)
235               if(imp_fcg) then
236           huaft(k)=coef1(k)*hu_mes(jm(k))+coef2(k)*hu_mes(jm(k)+1)
237           hvaft(k)=coef1(k)*hv_mes(jm(k))+coef2(k)*hv_mes(jm(k)+1)
238               endif ! imp_fcg
239               if(Turb_fcg) then
240           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))
241     $               +coef2(k)*hThTurb_mes(jm(k)+1)
242           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))
243     $               +coef2(k)*hqTurb_mes(jm(k)+1)
244               endif ! Turb_fcg
245             endif ! JM(k) .eq. 0
246            enddo
247            tsaft = ts_subr
248            pasprev=pas
249         else ! pas.gt.pasprev
250            print*,'timebef est:',timebef
251         endif ! pas.gt.pasprev    fin du bloc relatif au passage au pas
252                                  !de temps (meso) suivant
253c*** si on atteint le pas max des donnees experimentales ,on     ***
254c*** on conserve les derniers champs calcules                    ***
255      if(temps.ge.pasmax)then
256          do ll=1,klev
257               ht(ll)=htaft(ll)
258               hq(ll)=hqaft(ll)
259               hw(ll)=hwaft(ll)
260               hu(ll)=huaft(ll)
261               hv(ll)=hvaft(ll)
262               hThTurb(ll)=hThTurbaft(ll)
263               hqTurb(ll)=hqTurbaft(ll)
264          enddo
265               ts_subr = tsaft
266      else ! temps.ge.pasmax
267c*** on interpole sur les pas de temps de 10mn du gcm a partir   ***
268c** des pas de temps de 1h du meso_NH                            ***
269         do j=1,klev
270         ht(j)=((timeaft-time)*htbef(j)+(time-timebef)*htaft(j))/dt
271         hq(j)=((timeaft-time)*hqbef(j)+(time-timebef)*hqaft(j))/dt
272         hw(j)=((timeaft-time)*hwbef(j)+(time-timebef)*hwaft(j))/dt
273             if(imp_fcg) then
274         hu(j)=((timeaft-time)*hubef(j)+(time-timebef)*huaft(j))/dt
275         hv(j)=((timeaft-time)*hvbef(j)+(time-timebef)*hvaft(j))/dt
276             endif ! imp_fcg
277             if(Turb_fcg) then
278         hThTurb(j)=((timeaft-time)*hThTurbbef(j)
279     $           +(time-timebef)*hThTurbaft(j))/dt
280         hqTurb(j)= ((timeaft-time)*hqTurbbef(j)
281     $           +(time-timebef)*hqTurbaft(j))/dt
282             endif ! Turb_fcg
283         enddo
284              ts_subr = ((timeaft-time)*tsbef + (time-timebef)*tsaft)/dt
285       endif ! temps.ge.pasmax
286c
287        print *,' time,timebef,timeaft',time,timebef,timeaft
288        print *,' ht,htbef,htaft,hthturb,hthturbbef,hthturbaft'
289        do j= 1,klev
290           print *, j,ht(j),htbef(j),htaft(j),
291     $             hthturb(j),hthturbbef(j),hthturbaft(j)
292        enddo
293        print *,' hq,hqbef,hqaft,hqturb,hqturbbef,hqturbaft'
294        do j= 1,klev
295           print *, j,hq(j),hqbef(j),hqaft(j),
296     $             hqturb(j),hqturbbef(j),hqturbaft(j)
297        enddo
298c
299c-------------------------------------------------------------------
300c
301         IF (Ts_fcg) Ts = Ts_subr
302         return
303c
304c-----------------------------------------------------------------------
305c on sort les champs de "convergence" pour l instant initial 'in'
306c ceci se passe au pas temps itap=0 de la physique
307c-----------------------------------------------------------------------
308        entry get_uvd2(itap,dtime,file_forctl,file_fordat,
309     .           ht,hq,hw,hu,hv,hThTurb,hqTurb,ts,
310     .           imp_fcg,ts_fcg,Tp_fcg,Turb_fcg)
311             print*,'le pas itap est:',itap
312c
313c===================================================================
314c
315       write(*,'(a)') 'OPEN '//file_forctl
316       open(97,FILE=file_forctl,FORM='FORMATTED')
317c
318c------------------
319      do i=1,1000
320      read(97,1000,end=999) string
321 1000 format (a4)
322      if (string .eq. 'TDEF') go to 50
323      enddo
324 50   backspace(97)
325c-------------------------------------------------------------------
326c   *** on lit le pas de temps dans le fichier de donnees ***
327c   *** "forcing.ctl" et pasmax                           ***
328c-------------------------------------------------------------------
329      read(97,2000) aaa
330 2000  format (a80)
331         print*,'aaa est',aaa
332      aaa=spaces(aaa,1)
333         print*,'aaa',aaa
334      call getsch(aaa,' ',' ',5,atemps,nch)
335         print*,'atemps est',atemps
336        atemps=atemps(1:nch-2)
337         print*,'atemps',atemps
338        read(atemps,*) imn
339        dt=imn*60
340         print*,'le pas de temps dt',dt
341      call getsch(aaa,' ',' ',2,apasmax,nch)
342        apasmax=apasmax(1:nch)
343        read(apasmax,*) ipa
344        pasmax=ipa
345         print*,'pasmax est',pasmax
346      CLOSE(97)
347c------------------------------------------------------------------
348c *** on lit le pas de temps initial de la simulation ***
349c------------------------------------------------------------------
350                  in=itap
351cc                  pasprev=in
352cc                  time0=dt*(pasprev-1)
353                  pasprev=in-1
354                  time0=dt*pasprev
355C
356          close(98)
357c
358      write(*,'(a)') 'OPEN '//file_fordat
359      open(99,FILE=file_fordat,FORM='UNFORMATTED',
360     .          ACCESS='DIRECT',RECL=8)
361          icomp1 = nblvlm*4
362          IF (ts_fcg) icomp1 = icomp1 + 1
363          IF (imp_fcg) icomp1 = icomp1 + nblvlm*2
364          IF (Turb_fcg) icomp1 = icomp1 + nblvlm*2
365          icompt = icomp1*(in-1)
366          call rdgrads(99,icompt,nblvlm,z,ht_mes,hq_mes,hw_mes
367     .                   ,hu_mes,hv_mes,hthturb_mes,hqturb_mes
368     .                   ,ts_fcg,ts_subr,imp_fcg,Turb_fcg)
369          print *, 'get_uvd : rdgrads ->'
370          print *, tp_fcg
371c
372c following commented out because we have temperature already in ARM case
373c   (otherwise this is the potential temperature )
374c------------------------------------------------------------------------
375               if(Tp_fcg) then
376          do i = 1,nblvlm
377            ht_mes(i) = ht_mes(i)*(hplaym(i)/1000.)**rkappa
378          enddo
379               endif ! Tp_fcg
380           print*,'ht_mes ',(ht_mes(i),i=1,nblvlm)
381           print*,'hq_mes ',(hq_mes(i),i=1,nblvlm)
382           print*,'hw_mes ',(hw_mes(i),i=1,nblvlm)
383              if(imp_fcg) then
384           print*,'hu_mes ',(hu_mes(i),i=1,nblvlm)
385           print*,'hv_mes ',(hv_mes(i),i=1,nblvlm)
386           print*,'t',ts_subr
387              endif ! imp_fcg
388              if(Turb_fcg) then
389                 print*,'hThTurb_mes ',(hThTurb_mes(i),i=1,nblvlm)
390                 print*,'hqTurb ',     (hqTurb_mes(i),i=1,nblvlm)
391              endif ! Turb_fcg
392c----------------------------------------------------------------------
393c on a obtenu des champs initiaux sur les niveaux du meso_NH
394c on interpole sur les niveaux du gcm(niveau pression bien sur!)
395c-----------------------------------------------------------------------
396            do k=1,klev
397             if (JM(k) .eq. 0) then
398cFKC bug? ne faut il pas convertir tsol en tendance ????
399cRT bug taken care of by removing the stuff
400           htaft(k)=              ht_mes(jm(k)+1)
401           hqaft(k)=              hq_mes(jm(k)+1)
402           hwaft(k)=              hw_mes(jm(k)+1)
403               if(imp_fcg) then
404           huaft(k)=              hu_mes(jm(k)+1)
405           hvaft(k)=              hv_mes(jm(k)+1)
406               endif ! imp_fcg
407               if(Turb_fcg) then
408           hThTurbaft(k)=         hThTurb_mes(jm(k)+1)
409           hqTurbaft(k)=          hqTurb_mes(jm(k)+1)
410               endif ! Turb_fcg
411             else ! JM(k) .eq. 0
412           htaft(k)=coef1(k)*ht_mes(jm(k))+coef2(k)*ht_mes(jm(k)+1)
413           hqaft(k)=coef1(k)*hq_mes(jm(k))+coef2(k)*hq_mes(jm(k)+1)
414           hwaft(k)=coef1(k)*hw_mes(jm(k))+coef2(k)*hw_mes(jm(k)+1)
415               if(imp_fcg) then
416           huaft(k)=coef1(k)*hu_mes(jm(k))+coef2(k)*hu_mes(jm(k)+1)
417           hvaft(k)=coef1(k)*hv_mes(jm(k))+coef2(k)*hv_mes(jm(k)+1)
418               endif ! imp_fcg
419               if(Turb_fcg) then
420           hThTurbaft(k)=coef1(k)*hThTurb_mes(jm(k))
421     $                  +coef2(k)*hThTurb_mes(jm(k)+1)
422           hqTurbaft(k) =coef1(k)*hqTurb_mes(jm(k))
423     $                  +coef2(k)*hqTurb_mes(jm(k)+1)
424               endif ! Turb_fcg
425             endif ! JM(k) .eq. 0
426            enddo
427            tsaft = ts_subr
428c valeurs initiales des champs de convergence
429          do k=1,klev
430             ht(k)=htaft(k)
431             hq(k)=hqaft(k)
432             hw(k)=hwaft(k)
433                if(imp_fcg) then
434             hu(k)=huaft(k)
435             hv(k)=hvaft(k)
436                endif ! imp_fcg
437                if(Turb_fcg) then
438             hThTurb(k)=hThTurbaft(k)
439             hqTurb(k) =hqTurbaft(k)
440                endif ! Turb_fcg
441          enddo
442               ts_subr = tsaft
443          close(99)
444          close(98)
445c
446c-------------------------------------------------------------------
447c
448c
449 100      IF (Ts_fcg) Ts = Ts_subr
450        return
451c
452999     continue
453        stop 'erreur lecture, file forcing.ctl'
454        end
455
456      SUBROUTINE advect_tvl(dtime,zt,zq,vu_f,vv_f,t_f,q_f
457     :                     ,d_t_adv,d_q_adv)
458      use dimphy
459      implicit none
460
461#include "dimensions.h"
462ccccc#include "dimphy.h"
463
464      integer k
465      real dtime, fact, du, dv, cx, cy, alx, aly
466      real zt(klev), zq(klev,3)
467     :   , vu_f(klev), vv_f(klev), t_f(klev), q_f(klev,3)
468
469      real d_t_adv(klev), d_q_adv(klev,3)
470
471c Velocity of moving cell
472      data cx,cy /12., -2./
473
474c Dimensions of moving cell
475      data alx,aly /100 000.,150 000./
476
477      do k = 1, klev
478            du = abs(vu_f(k)-cx)/alx
479            dv = abs(vv_f(k)-cy)/aly
480            fact = dtime *(du+dv-du*dv*dtime)
481            d_t_adv(k) = fact * (t_f(k)-zt(k))
482            d_q_adv(k,1) = fact * (q_f(k,1)-zq(k,1))
483            d_q_adv(k,2) = fact * (q_f(k,2)-zq(k,2))
484            d_q_adv(k,3) = fact * (q_f(k,3)-zq(k,3))
485      enddo
486
487      return
488      end
489
490      SUBROUTINE copie(klevgcm,playgcm,psolgcm,file_forctl)
491      implicit none
492
493ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
494c cette routine remplit les COMMON com1_phys_gcss et com2_phys_gcss.h
495cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
496
497      INTEGER klev !nombre de niveau de pression du GCM
498      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
499      INTEGER JM(100)
500      REAL coef1(100)   !coefficient d interpolation
501      REAL coef2(100)   !coefficient d interpolation
502
503      INTEGER nblvlm !nombre de niveau de pression du mesoNH
504      REAL playm(100) !pression en Pa au milieu de chaque couche Meso-NH
505      REAL hplaym(100)!pression en hecto-Pa des milieux de couche Meso-NH
506
507      COMMON/com1_phys_gcss/klev,play,JM,coef1,coef2
508      COMMON/com2_phys_gcss/nblvlm,playm,hplaym
509
510      integer i,k,klevgcm
511      real playgcm(klevgcm) ! pression en milieu de couche du gcm
512      real psolgcm
513      character*80 file_forctl
514
515      klev = klevgcm
516
517c---------------------------------------------------------------------
518c pression au milieu des couches du gcm dans la physiq
519c (SB: remplace le call conv_lipress_gcm(playgcm) )
520c---------------------------------------------------------------------
521
522       do k = 1, klev
523        play(k) = playgcm(k)
524        print*,'la pression gcm est:',play(k)
525       enddo
526
527c----------------------------------------------------------------------
528c lecture du descripteur des donnees Meso-NH (forcing.ctl):
529-> nb niveaux du meso.NH (nblvlm) + pressions meso.NH
530c (on remplit le COMMON com2_phys_gcss)
531c----------------------------------------------------------------------
532
533      call mesolupbis(file_forctl)
534
535      print*,'la valeur de nblvlm est:',nblvlm
536
537c----------------------------------------------------------------------
538c etude de la correspondance entre les niveaux meso.NH et GCM;
539c calcul des coefficients d interpolation coef1 et coef2
540c (on remplit le COMMON com1_phys_gcss)
541c----------------------------------------------------------------------
542
543      call corresbis(psolgcm)
544
545c---------------------------------------------------------
546c TEST sur le remplissage de com1_phys_gcss et com2_phys_gcss:
547c---------------------------------------------------------
548 
549      write(*,*) ' '
550      write(*,*) 'TESTS com1_phys_gcss et com2_phys_gcss dans copie.F'
551      write(*,*) '--------------------------------------'
552      write(*,*) 'GCM: nb niveaux:',klev,' et pression, coeffs:'
553      do k = 1, klev
554      write(*,*) play(k), coef1(k), coef2(k)
555      enddo
556      write(*,*) 'MESO-NH: nb niveaux:',nblvlm,' et pression:'
557      do k = 1, nblvlm
558      write(*,*) playm(k), hplaym(k)
559      enddo
560      write(*,*) ' '
561
562      end
563      SUBROUTINE mesolupbis(file_forctl)
564      implicit none
565c
566cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
567c
568c Lecture descripteur des donnees MESO-NH (forcing.ctl):
569c -------------------------------------------------------
570c
571c     Cette subroutine lit dans le fichier de controle "essai.ctl"
572c     et affiche le nombre de niveaux du Meso-NH ainsi que les valeurs
573c     des pressions en milieu de couche du Meso-NH (en Pa puis en hPa).
574cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
575c
576      INTEGER nblvlm !nombre de niveau de pression du mesoNH
577      REAL playm(100)  !pression en Pa milieu de chaque couche Meso-NH
578      REAL hplaym(100) !pression en hPa des milieux de couche Meso-NH
579      COMMON/com2_phys_gcss/nblvlm,playm,hplaym
580
581      INTEGER i,lu,k,mlz,mlzh,j
582 
583      character*80 file_forctl
584
585      character*4 a
586      character*80 aaa,anblvl,spaces
587      integer nch
588
589      lu=9
590      open(lu,file=file_forctl,form='formatted')
591c
592      do i=1,1000
593      read(lu,1000,end=999) a
594      if (a .eq. 'ZDEF') go to 100
595      enddo
596c
597 100  backspace(lu)
598      print*,'  DESCRIPTION DES 2 MODELES : '
599      print*,' '
600c
601      read(lu,2000) aaa
602 2000  format (a80)
603       aaa=spaces(aaa,1)
604       call getsch(aaa,' ',' ',2,anblvl,nch)
605         read(anblvl,*) nblvlm
606
607c
608      print*,'nbre de niveaux de pression Meso-NH :',nblvlm
609      print*,' '
610      print*,'pression en Pa de chaque couche du meso-NH :'
611c
612      read(lu,*) (playm(mlz),mlz=1,nblvlm)
613c      Si la pression est en HPa, la multiplier par 100
614      if (playm(1) .lt. 10000.) then
615        do mlz = 1,nblvlm
616         playm(mlz) = playm(mlz)*100.
617        enddo
618      endif
619      print*,(playm(mlz),mlz=1,nblvlm)
620c
621 1000 format (a4)
622 1001 format(5x,i2)
623c
624      print*,' '
625      do mlzh=1,nblvlm
626      hplaym(mlzh)=playm(mlzh)/100.
627      enddo
628c
629      print*,'pression en hPa de chaque couche du meso-NH: '
630      print*,(hplaym(mlzh),mlzh=1,nblvlm)
631c
632      close (lu)
633      return
634c
635 999  stop 'erreur lecture des niveaux pression des donnees'
636      end
637
638      SUBROUTINE rdgrads(itape,icount,nl,z,ht,hq,hw,hu,hv,hthtur,hqtur,
639     :  ts_fcg,ts,imp_fcg,Turb_fcg)
640      IMPLICIT none
641      INTEGER itape,icount,icomp, nl
642      real z(nl),ht(nl),hq(nl),hw(nl),hu(nl),hv(nl)
643      real hthtur(nl),hqtur(nl)
644      real ts
645c
646      INTEGER i, k
647c
648      LOGICAL imp_fcg,ts_fcg,Turb_fcg
649c
650      icomp = icount
651c
652c
653         do k=1,nl
654            icomp=icomp+1
655            read(itape,rec=icomp)z(k)
656            print *,'icomp,k,z(k) ',icomp,k,z(k)
657         enddo
658         do k=1,nl
659            icomp=icomp+1
660            read(itape,rec=icomp)hT(k)
661             print*, hT(k), k
662         enddo
663         do k=1,nl
664            icomp=icomp+1
665            read(itape,rec=icomp)hQ(k)
666         enddo
667c
668             if(turb_fcg) then
669         do k=1,nl
670            icomp=icomp+1
671           read(itape,rec=icomp)hThTur(k)
672         enddo
673         do k=1,nl
674            icomp=icomp+1
675           read(itape,rec=icomp)hqTur(k)
676         enddo
677             endif
678         print *,' apres lecture hthtur, hqtur'
679c
680          if(imp_fcg) then
681
682         do k=1,nl
683            icomp=icomp+1
684           read(itape,rec=icomp)hu(k)
685         enddo
686         do k=1,nl
687            icomp=icomp+1
688            read(itape,rec=icomp)hv(k)
689         enddo
690
691          endif
692c
693         do k=1,nl
694            icomp=icomp+1
695            read(itape,rec=icomp)hw(k)
696         enddo
697c
698              if(ts_fcg) then
699         icomp=icomp+1
700         read(itape,rec=icomp)ts
701              endif
702c
703      print *,' rdgrads ->'
704
705      RETURN
706      END
707
708      SUBROUTINE corresbis(psol)
709      implicit none
710
711ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
712c Cette subroutine calcule et affiche les valeurs des coefficients
713c d interpolation qui serviront dans la formule d interpolation elle-
714c meme.
715cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
716
717      INTEGER klev    !nombre de niveau de pression du GCM
718      REAL play(100)  !pression en Pa au milieu de chaque couche GCM
719      INTEGER JM(100)
720      REAL coef1(100) !coefficient d interpolation
721      REAL coef2(100) !coefficient d interpolation
722
723      INTEGER nblvlm !nombre de niveau de pression du mesoNH
724      REAL playm(100) !pression en Pa milieu de chaque couche Meso-NH
725      REAL hplaym(100)!pression en hPa des milieux de couche Meso-NH
726
727      COMMON/com1_phys_gcss/klev,play,JM,coef1,coef2
728      COMMON/com2_phys_gcss/nblvlm,playm,hplaym
729
730      REAL psol
731      REAL val
732      INTEGER k, mlz, mlzh
733
734
735      do k=1,klev
736         val=play(k)
737       if (val .gt. playm(1)) then
738          mlz = 0
739          JM(1) = mlz
740          coef1(1)=(playm(mlz+1)-val)
741     *             /(playm(mlz+1)-psol)
742          coef2(1)=(val-psol)
743     *             /(playm(mlz+1)-psol)
744       else if (val .gt. playm(nblvlm)) then
745         do mlz=1,nblvlm
746          if (     val .le. playm(mlz)
747     *       .and. val .gt. playm(mlz+1))then
748           JM(k)=mlz
749           coef1(k)=(playm(mlz+1)-val)
750     *              /(playm(mlz+1)-playm(mlz))
751           coef2(k)=(val-playm(mlz))
752     *              /(playm(mlz+1)-playm(mlz))
753          endif
754         enddo
755       else
756         JM(k) = nblvlm-1
757         coef1(k) = 0.
758         coef2(k) = 0.
759       endif
760      enddo
761c
762cc      if (play(klev) .le. playm(nblvlm)) then
763cc         mlz=nblvlm-1
764cc         JM(klev)=mlz
765cc         coef1(klev)=(playm(mlz+1)-val)
766cc     *            /(playm(mlz+1)-playm(mlz))
767cc         coef2(klev)=(val-playm(mlz))
768cc     *            /(playm(mlz+1)-playm(mlz))
769cc      endif
770c
771      print*,' '
772      print*,'         INTERPOLATION  : '
773      print*,' '
774      print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
775      print*,(JM(k),k=1,klev)
776      print*,'correspondance de 9 niveaux du GCM sur les 53 du meso-NH:'
777      print*,(JM(k),k=1,klev)
778      print*,' '
779      print*,'valeurs du premier coef d"interpolation pour les 9 niveaux
780     *: '
781      print*,(coef1(k),k=1,klev)
782      print*,' '
783      print*,'valeurs du deuxieme coef d"interpolation pour les 9 niveau
784     *x: '
785      print*,(coef2(k),k=1,klev)
786c
787      return
788      end
789      SUBROUTINE GETSCH(STR,DEL,TRM,NTH,SST,NCH)
790C***************************************************************
791C*                                                             *
792C*                                                             *
793C* GETSCH                                                      *
794C*                                                             *
795C*                                                             *
796C* modified by :                                               *
797C***************************************************************
798C*   Return in SST the character string found between the NTH-1 and NTH
799C*   occurence of the delimiter 'DEL' but before the terminator 'TRM' in
800C*   the input string 'STR'. If TRM=DEL then STR is considered unlimited.
801C*   NCH=Length of the string returned in SST or =-1 if NTH is <1 or if
802C*   NTH is greater than the number of delimiters in STR.
803      IMPLICIT INTEGER (A-Z)
804      CHARACTER STR*(*),DEL*1,TRM*1,SST*(*)
805      NCH=-1
806      SST=' '
807      IF(NTH.GT.0) THEN
808        IF(TRM.EQ.DEL) THEN
809          LENGTH=LEN(STR)
810        ELSE
811          LENGTH=INDEX(STR,TRM)-1
812          IF(LENGTH.LT.0) LENGTH=LEN(STR)
813        ENDIF
814C*     Find beginning and end of the NTH DEL-limited substring in STR
815        END=-1
816        DO 1,N=1,NTH
817        IF(END.EQ.LENGTH) RETURN
818        BEG=END+2
819        END=BEG+INDEX(STR(BEG:LENGTH),DEL)-2
820        IF(END.EQ.BEG-2) END=LENGTH
821C*        PRINT *,'NTH,LENGTH,N,BEG,END=',NTH,LENGTH,N,BEG,END
822    1   CONTINUE
823        NCH=END-BEG+1
824        IF(NCH.GT.0) SST=STR(BEG:END)
825      ENDIF
826      END
827      CHARACTER*(*) FUNCTION SPACES(STR,NSPACE)
828C
829C CERN PROGLIB# M433    SPACES          .VERSION KERNFOR  4.14  860211
830C ORIG.  6/05/86 M.GOOSSENS/DD
831C
832C-    The function value SPACES returns the character string STR with
833C-    leading blanks removed and each occurence of one or more blanks
834C-    replaced by NSPACE blanks inside the string STR
835C
836      CHARACTER*(*) STR
837C
838      LENSPA = LEN(SPACES)
839      SPACES = ' '
840      IF (NSPACE.LT.0) NSPACE = 0
841      IBLANK = 1
842      ISPACE = 1
843  100 INONBL = INDEXC(STR(IBLANK:),' ')
844      IF (INONBL.EQ.0) THEN
845          SPACES(ISPACE:) = STR(IBLANK:)
846                                                    GO TO 999
847      ENDIF
848      INONBL = INONBL + IBLANK - 1
849      IBLANK = INDEX(STR(INONBL:),' ')
850      IF (IBLANK.EQ.0) THEN
851          SPACES(ISPACE:) = STR(INONBL:)
852                                                    GO TO 999
853      ENDIF
854      IBLANK = IBLANK + INONBL - 1
855      SPACES(ISPACE:) = STR(INONBL:IBLANK-1)
856      ISPACE = ISPACE + IBLANK - INONBL + NSPACE
857      IF (ISPACE.LE.LENSPA)                         GO TO 100
858  999 END
859      FUNCTION INDEXC(STR,SSTR)
860C
861C CERN PROGLIB# M433    INDEXC          .VERSION KERNFOR  4.14  860211
862C ORIG. 26/03/86 M.GOOSSENS/DD
863C
864C-    Find the leftmost position where substring SSTR does not match
865C-    string STR scanning forward
866C
867      CHARACTER*(*) STR,SSTR
868C
869      LENS   = LEN(STR)
870      LENSS  = LEN(SSTR)
871C
872      DO 10 I=1,LENS-LENSS+1
873          IF (STR(I:I+LENSS-1).NE.SSTR) THEN
874              INDEXC = I
875                                         GO TO 999
876          ENDIF
877   10 CONTINUE
878      INDEXC = 0
879C
880  999 END
Note: See TracBrowser for help on using the repository browser.