source: LMDZ5/branches/testing/libf/phylmd/dyn1d/1Dconv.h @ 2408

Last change on this file since 2408 was 2408, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2298:2396 into testing branch

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